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Chapter 1 
Introduction 


This is the final report for the NASA funded project NAG-1-1184 entitled “Crack 
Growth Prediction Methodology for Multi-Site Damage.” The primary objective of 
the project was to create a capability to simulate curvilinear fatigue crack growth 
and ductile tearing in aircraft fuselages subjected to widespread fatigue damage. 
The second objective was to validate the capability by way of comparisons to ex- 
perimental results. Both objectives have been achieved and the results are detailed 
herein. The body of this report is derived primarily from the Ph.D. thesis of the 
first author [19]. 

1.1 Overview 

Modern aircraft structures are designed using a damage tolerance philosophy. This 
design philosophy envisions sufficient strength and structural integrity of the air- 
craft to sustain major damage and to avoid catastrophic failure. However, struc- 
tural aging of the aircraft may significantly reduce the strength below an acceptable 
level; this raises many important safety issues. 

Concerns about aging aircraft are reinforced by the in-flight structural failure 
of an Aloha Airlines Boeing 737 on April 28, 1988 [97]. The failure precipitated 
from the link-up of small fatigue cracks extending from adjacent rivet holes in a 
fuselage lap-splice joint. This caused approximately 18 feet of the upper crown 
skin and structure to separate from the fuselage (see Figures 1.1 and 1.2). The 
1988 Aloha Airlines accident created a revolution in the aircraft community. The 
problems associated with aging aircraft have to be quantified and the methodology 
to ensure the structural integrity of airplanes has to be reassessed [5]. 

One of the most notable problems in aging aircraft is widespread fatigue damage 
(WFD) defined in [139] as “the simultaneous presence of fatigue cracks at multi- 
ple structural details that are of sufficient size and density whereby the structure 
will no longer meet its damage tolerance requirement.” In response to such aging 
aircraft problems, the National Aeronautics and Space Administration (NASA) ini- 
tiated an Airframe Structural Integrity Program (ASIP) [52, 95, 53, 54] to develop 
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Figure 1.1: 1988 Aloha Airlines accident (aircraft after landing). 



Figure 1.2: 1988 Aloha Airlines accident-the shaded area illustrates the part of 
fuselage lost at cruising altitude. 
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an advanced analysis methodology for predicting structural integrity and residual 
strength 1 of fuselage structures with WFD. The analysis methodology would al- 
low engineers to maintain the aging fleet economically while insuring continuous 
airworthiness. 

The work described in this thesis is part of the ASIP program. The main 
objective of this work is to develop a computational analysis methodology to sim- 
ulate realistic crack growth and to predict remaining life and residual strength of 
complex built-up structures. While the methodology developed is generic in na- 
ture, the particular focus is on fuselage structural integrity where WFD is likely 
to occur as the fleet grows older. The analysis methodology will help to determine 
service inspection intervals, quantitatively evaluate inspection findings, and design 
and certify damage-tolerant structural repairs. Thus, the outcome will improve 
the technology to support the safe operation of the current fleet and the design of 
more damage-tolerant aircraft for the next-generation fleet. 

1.2 Background 1: Structural Integrity of Fuse- 
lage Structures 

A general overview of the structural integrity of a pressurized fuselage with cracks 
is provided as background material for the thesis. In particular, the following issues 
are discussed: 

• concerns about WFD related to the loss of residual strength, 

• the characterization of WFD in riveted fuselage structures, 

• typical life of aircraft structures and the dominant behavior of cracks, and 
finally 

• the analysis methodology described in this thesis to evaluate the structural 
integrity of damaged structures. 

1.2.1 WFD in Fuselage Structures 

The philosophy of damage tolerance presumes that any damage initiated by fatigue, 
accident, or corrosion will be found before catastrophic failure [138]. The safety 
of the aircraft heavily depends upon finding cracks before they reach a critical 
length. The occurrence of WFD, however, drastically reduces the residual strength 
or decreases the critical crack size as illustrated in Figure 1.3. The loss of residual 
strength in the presence of WFD has raised great safety concerns for aging aircraft. 


1 Residual strength is the maximum load carrying capacity of a damaged struc- 
ture. 
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Figure 1.3: Illustration of effect of WFD on residual strength and critical crack 
size. 


To establish the characteristics of WFD in fuselage riveted structures, Piascik 
et al. [105] and Harris et al. [54] conducted teardown fractographic examinations 
of aircraft components. A four bay section of a Boeing 747 fuselage containing 
a longitudinal lap splice joint was examined after conducting a full-scale fatigue 
test with 60,000 full pressurization cycles; this is about three times the original 
economic design life of the aircraft. Several observations were made from detailed 
non-destructive and destructive examinations of each rivet hole: 

1. crack initiation mechanisms included high local stresses, fretting along mat- 
ing surfaces, and manufacturing defects created during the riveting process; 

2. fatigue cracks were present at virtually every rivet hole in the upper row of 
the lap joint; 

3. the lengths of all of the fatigue cracks at link-up were approximately the 
same. 

The last observation implies that as long as the crack has extended a consider- 
able distance from the rivet head, the crack growth behavior is somewhat indepen- 
dent of the initiating mechanism. Thus, the typical life of aircraft structures can 
be subdivided into the nucleation and crack growth periods as shown in Figure 1.4. 
The nucleation period heavily depends on micro-structural details of the material. 
The microscopic studies provide fundamental understanding and phenomenological 
criteria for fatigue and fracture used in macro-scale applications. Regardless of the 
initiation mechanism, fracture mechanics is adequate to describe the macrocrack 
behavior for practical problems. The methodology developed herein is intended 
for the macro-scale and it is the crack growth period that is of primary interest. 
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nucleation period crack growth period 

Figure 1.4: Illustration of typical life of aircraft structures. 

1.2.2 Crack Growth in Thin Sheet Metals 

In general, crack growth in thin, ductile materials used in aircraft fuselages is likely 
to experience fatigue and stable crack growth before the occurrence of fast frac- 
ture and final failure. Fatigue is a process of cumulative damage that is caused 
by repeated fluctuating loads. Fatigue crack propagation can be characterized by 
a crack growth-rate model that predicts the number of loading cycles required to 
propagate a fatigue crack to a critical size. Stress intensity factors (SIFs) un- 
der fatigue loading are below the critical value for quasi-static or unstable crack 
propagation. Under such circumstance, linear elastic fracture mechanics (LEFM) 
suffices to characterize the crack growth-rate model. Stable crack growth and final 
failure generally occur at the very last loading cycle of the life of aircraft. Crack 
propagation at this stage involves elastic-plastic stable tearing followed by fast- 
fracture. Since crack growth is no longer under small-scale yielding conditions, 
elastic-plastic fracture mechanics (EPFM) is needed to characterize the fracture 
behavior and to predict the residual strength. 

1.2.3 Bulging in Pressurized Thin Shell Structures 

For cracks in a pressurized fuselage, the out-of-plane deformation or bulging at 
crack edges is an essential characteristic feature of the displacement fields. Fol- 
lowing Potyondy et al. [108], the dominant factors that affect the behavior of 
through-cracks in the skin of pressurized fuselages are: 

1. a geometrically nonlinear stiffening effect that restricts the crack edge bulging, 

2. the presence of stiffening elements that alter the stress distribution in the 
skin, 

3. the internal pressure and the mechanical loads that act on the structure, and 

4. plasticity effects. 

All the factors described above are taken into consideration in developing the 
analysis methodology in this work. The following general guidelines are used to 
characterize crack growth and to evaluate fuselage structural integrity: 
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• For fatigue crack growth, use thin shell analyses with geometric nonlinearity 
to evaluate SIFs for shells under membrane and bending loading. The SIFs 
are then used to evaluate crack growth-rate. 

• For stable crack growth, use thin shell analyses with geometric and material 
nonlinearity to evaluate crack tip opening angle (CTOA). The critical CTOA 
is then used to control the crack advancement. Structural integrity and 
residual strength are then evaluated based on the predicted results of elastic- 
plastic crack growth simulations. 


1.3 Background 2: Computational Issues for Crack 
Growth Simulation 

To simulate realistic crack growth where crack trajectories are not known a priori , 
continual updating of the geometry is required. This feature makes conventional 
programs for computational solid mechanics difficult to use, if used alone. In this 
section, a brief overview of the main computational issues for arbitrary, discrete 
crack growth simulations is provided. Key aspects of the high-level description of 
the implementation to make the crack growth simulations efficient are presented 2 . 

Discrete crack growth simulation is an incremental process, where a series of 
steps is repeated for a progression of models. The process continues until a suitable 
termination condition is reached. Results of such a simulation might include one or 
more of the following: a final crack geometry, a loading versus crack size history, 
a crack opening profile, or a history of the crack-front, fracture parameters. In 
general, each increment in the process relies on previously computed results and 
represents one crack configuration. Following Carter et al. [17, 18], data in each 
crack growth increment can be divided into: representation database (R;, where 
the subscript denotes the increment number), analysis database (Ai), equilibrium 
database (Ei), and fracture parameter database (F;). Each simulation of crack 
growth increment involves three major processes, a discretization process (D), a 
solution process (S), and an update process (U). The sequence of operations is 
illustrated in Figure 1.5. A discretization process, D, primarily consisting of a 
meshing function, M, transforms a representational description of a cracked body 
to a discrete model suitable for stress analysis. A solution process, S, computes 
unknown field variables, E ; , and fracture parameters, F;. An update process, U, 
takes the equilibrium state field variables and the existing representation, using 
a function that predicts crack shape evolution, C, creates a new representational 
database. The major problem in using conventional programs alone to simulate 

2 It is not the intent of the author to describe details of computational issues 
about arbitrary crack growth, but instead present a broad background and a de- 
parture point to introduce the software written and developed as part of this thesis. 
The readers are referred to [143, 17, 18] for more detailed and in-depth discussion. 
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Figure 1.5: Illustration of the sequence of operations on databases and processes 
involved in each crack growth increment. 


discrete crack growth is that they are ill-suited for tracking geometry changes that 
accompany crack growth unless the growth is constrained to be along existing 
element edges or faces. For general cases, where crack trajectories are not known 
a priori , a general representational model provides the key to efficiently handle 
crack growth simulations. 

1.4 FRANC3D/STAGS Software Environment 

A software system FRANC3D/STAGS that incorporates the conceptual crack 
growth model described above to support the analysis methodology in evaluating 
structural integrity of aging aircraft is developed. The system combines a topology 
based program, FRANC3D (FRacture ANalysis Code for 3D solids and shells), and 
a general nonlinear shell finite element program, STAGS (STructual Analysis of 
General Shells). FRANC3D has been under development by the Cornell Fracture 
Group since the late 1980s. The aim is a systematic representational model for 
arbitrary crack growth analysis. The key concepts embodied in FRANC3D are: 

• solid modeling tools, 

• a topological data structure that separates topology from geometry, 

• the association of model attributes with topological entities, and 
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• a hierarchy of topological models to organize and guide the discretization 
process. 

The representational model in FRANC3D is well-suited for tracking the geom- 
etry changes that accompany crack growth. The task of updating the representa- 
tional model and generating a sequence of analysis models is separated from the 
task of performing an analysis on a particular model. For the FRANC3D/STAGS 
software, the analyses are performed by the STAGS code. 

STAGS is a finite element code for general-purpose analysis of shell structures 
developed by Lockheed-Martin’s Advanced Technology Center [113]. The main 
analysis capabilities in STAGS are: 

• linear elastic stress analysis, 

• geometrically and materially nonlinear stress analysis, 

• linear bifurcation buckling analysis, and 

• transient response analysis. 

The FRANC3D/STAGS software first developed by Potyondy [106] is further 
modified by the author to support the evolving analysis methodology for evaluating 
structural integrity of aircraft structures. 

The software components of FRANC3D/STAGS necessary to perform arbitrary 
crack growth simulation are illustrated in Figure 1.6. The FRANC3D code con- 
trols the entire process, allowing the analyst to compute the equilibrium states for 
a series of structural configurations. The STAGS code performs the stress analy- 
sis. Finally, an interface program is written to facilitate the data communication 
between the two codes. 


1.5 Organization of the Dissertation 

The primary objective of the dissertation is to develop an accurate structural 
analysis methodology and a useful and usable software program for predicting the 
structural integrity and residual strength of fuselage structures. The dissertation 
is divided into two parts: (1) elastic-plastic crack growth analyses and residual 
strength prediction with self-similar crack growth, Chapters 2-4; and (2) crack 
trajectory prediction with non-self-similar, curvilinear crack growth, Chapters 5- 
7. 

Chapter 2 reviews and critiques various fracture mechanics methods for sim- 
ulating elastic-plastic crack growth and predicting residual strength of thin-sheet 
metallic structures. Among the methods, the CTOA fracture criterion is found 
to be superior due to its relative independence of the geometry of the structure, 
the length of the crack, and the presence of multiple cracks. The concepts and 
formulations of the CTOA criterion are presented. Elastic-plastic crack growth, 
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Figure 1.6: The FRANC3D/STAGS software components, their primary functions, 
and their interactions necessary to perform arbitrary crack growth sim- 
ulation (modified after [106]). 
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link-up of multiple cracks, and residual strength analyses using the CTOA fracture 
criterion are discussed. 

Chapters 3 and 4 describe finite element analyses of fracture tests using the 
CTOA fracture criterion. Elastic-plastic crack growth simulations and residual 
strength prediction of fracture coupon tests and full-scale fuselage panel tests are 
presented. In Chapter 3, fracture behavior in coupon tests including middle-crack 
tension (MT) and multiple crack specimens is analyzed. In Chapter 4, fracture 
behavior in full-scale narrow body and wide body panel tests is analyzed. Possible 
scenarios that can occur in pressurized fuselages are examined, including lead crack 
growth, multi-site damage, multiple crack interaction, plastic wake from fatigue 
crack growth, tear strap failure, and corrosion damage. 

Chapters 5 through 7, the second part of the dissertation, consist of mate- 
rial related to crack trajectory prediction. An evolving methodology to improve 
structural integrity of aircraft structures using the crack turning phenomenon is 
discussed. 

Chapter 5 reviews and critiques various crack growth directional criteria for 
crack trajectory prediction. A directional criterion based on the maximum tan- 
gential stress theory, but taking into account the effect of T-stress and fracture 
toughness orthotropy is developed. In Chapter 6, the path independent contour 
integral method for T-stress evaluation is presented. The numerical accuracy using 
the path independent integral is assessed by highly accurate two-dimensional p- 
and hp-version adaptive finite element analyses. Chapter 7 analyzes curvilinear 
crack growth in coupon tests and in full-scale narrow body fuselage panel tests. 
The T-stress and fracture toughness orthotropy effect on crack trajectory predic- 
tion is examined. 

The final chapter summarizes the contributions of this thesis, draws conclu- 
sions, and where appropriate, provides recommendations for future work. 



Chapter 2 

Theory for CTOA-Driven 
Elastic-Plastic Crack Growth and 
Residual Strength Analysis 


This chapter together with Chapters 3 and 4 gives a comprehensive treatise on us- 
ing the crack tip opening angle (CTOA) fracture criterion to predict elastic-plastic 
crack growth and residual strength of thin-sheet metallic structures. Theories, con- 
cepts, and formulations related to the CTOA-driven fracture criterion are given in 
this chapter. Elastic-plastic crack growth simulations and residual strength pre- 
diction of coupon tests and full-scale fuselage panel tests are presented in the next 
two chapters. 

2.1 Introduction 

To predict successfully fracture behavior and residual strength of aircraft fuselage 
structures subjected to widespread fatigue damage (WFD), a fracture criterion 
independent of the geometry of the structure, the length of the crack, and the 
presence of multi-site damage (MSD) is required [29]. In this chapter, a brief 
evaluation of various fracture mechanics methods to simulate elastic-plastic crack 
growth and to predict residual strength of damaged structures is given. The evalu- 
ation focuses on the applicability to thin-sheet metallic structures with single and 
multiple cracks where plastic flow makes a substantial contribution to crack growth 
resistance. 

Among the fracture methods, the superior nature of the CTOA fracture cri- 
terion to characterize elastic-plastic crack growth in thin-sheet metals is revealed 
after review of two recent evaluations [91, 28]. Theories for simulating the CTOA- 
driven crack growth are discussed and guidelines for using the CTOA fracture 
criterion to predict residual strength of aircraft structures are presented. 
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2.2 Evaluation of Fracture Mechanics Methods 

Thin sheet metallic structures in ductile states generally undergo stable crack 
growth before the occurrence of fast fracture. To have an accurate and reliable 
prediction of residual strength for such structures, a fracture criterion that can 
characterize stable crack growth under conditions of general yielding is needed. 
Various fracture mechanics methods have been developed over the past several 
decades. However, for materials exhibiting a large amount of plasticity and stable 
crack growth prior to failure, there is no consensus on the most satisfactory method 
[116, 91]. To evaluate various methods in assessing crack growth resistance and 
predicting failure of flawed structures, an experimental and predictive round robin 
was conducted in 1979-1980 by the ASTM Committee E-24 on Fracture Testing 
[91]. The fracture analysis methods used in the round robin included: 

1. linear elastic fracture mechanics (LEFM) corrected for size effects or for 
plastic yielding [91], 

2. equivalent energy [148], 

3. the two-parameter fracture criterion, I\' F and m [87], 

4. the deformation plasticity failure assessment diagram based on deformation 
plasticity, a ./-integral estimation scheme, and a solution from the Plastic 
Handbook [10, 75], 

5. the theory of ductile fracture [11], 

6. the K r - curve with the Dugdale model [37], 

7. the effective A#-curve [79], 

8. a two-dimensional (2D) finite element analysis using the CTOA criterion 
with stable crack growth [90], and 

9. a three-dimensional (3D) finite element analysis using a crack-front singular- 
ity parameter with a stationary crack [81]. 

Fracture tests were conducted on compact tension specimens (CT), middle- 
crack tension specimens (MT), and three-hole-crack tension specimens (THCT) 
as shown in Figure 2.1. Three materials tested were 7075-T651 aluminum alloy, 
2024-T351 aluminum alloy, and 304 stainless steel. The accuracy of the prediction 
methods was judged by the failure loads obtained from experiments. For 7075- 
T651 aluminum alloy, the best methods were the effective A'^-curve, a 2D finite 
element analysis using CTOA with stable crack growth, and the AO-curve with 
the Dugdale model. For 2024-T351 aluminum alloy, the best methods were the 
two-parameter fracture criterion, a 2D finite element analysis using CTOA with 
stable crack growth, the /Or curve with the Dugdale model, the effective /Or curve, 
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Figure 2.1: Specimen configurations tested and analyzed in ASTM round robin 
(after [91]). 

and the deformation plasticity failure assessment diagram. For 304 stainless steel, 
the best methods were the effective A'^-curve, a 2D finite element analysis using 
CTOA with stable crack growth, the two-parameter fracture criterion, and the 
deformation plasticity failure assessment diagram. 

These tests were conducted using specimens with a single crack configuration. 
Recently, fracture tests were conducted on a 0.09 inch thick, 2024-T3 aluminum 
alloy using CT, MT, and MSD specimens [28]. Several fracture mechanics methods 
were again evaluated including the effective K R - curve [79, 3], the J-integral resis- 
tance curve {J R ) [92], the crack-opening resistance curve {Sr) [55], the T*-integral 
resistance curve {T R ) [4], the plastic-zone link-up criterion [139], and the critical 
CTOA fracture criterion using a 3D finite element analysis [29]. The study con- 
cluded that the plastic-zone link-up criterion had limited use in predicting fracture 
behavior of specimens. The effective K R , J R , T^, and Sr fracture criteria could 
predict MSD fracture behavior of some larger specimens based on small specimen 
tests, but were limited to a certain size of specimens. The critical CTOA fracture 
criterion using a 3D, elastic-plastic finite element analysis was able to predict the 
fracture behavior for all specimen sizes. 

Based on the above evaluations, the CTOA-driven, elastic-plastic stable crack 
growth simulation appears to be a plausible fracture analysis method to assess 
crack growth resistance and to predict residual strength of thin-sheet metallic 
structures. 
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2.3 CTOA Fracture Criterion 

The CTOA fracture criterion is essentially an integration of the near-tip strains. It 
evolves from the critical crack tip opening displacement (CTOD) concept proposed 
by Wells [144], Since the CTOD has a limiting value of zero at the crack tip, the 
local slope of the crack tip opening profile, or CTOA, was suggested to characterize 
the crack growth behavior [2, 36]. Newman [90], Rice and Sorensen [121], and 
Kanninen and Popelar [64] further defined the CTOA as the crack tip opening 
angle measured at a fixed distance behind the moving crack tip. 

The CTOA fracture criterion asserts that the angle maintains a constant value 
during stable crack growth for a given thickness of a metallic material. This phe- 
nomenon has been observed in numerous experiments for a wide range of metals 
[66, 63, 65, 90, 34, 94], indirectly supported from slip-line field solutions [121, 118], 
and verified by numerical simulations [36, 90, 64, 93, 35, 30, 27, 29, 28, 20, 23, 22], 
Tests on aluminum alloys as well as steels [65, 90] have confirmed that the CTOA 
is essentially constant after a certain transitional period of stable crack growth. A 
larger critical CTOA during the initiation of stable tearing rapidly decreases to a 
constant value. The amount of crack growth to reach the constant CTOA is ap- 
proximately equal to the specimen thickness [34], Based on a fatigue marker load 
technique and scanning electron microscope observations, Dawicke and Sutton [34] 
concluded that the non-constant CTOA region is associated with severe tunneling 
during the initiation of stable crack growth. 

Asymptotic solutions of a growing crack provide indirect support for using the 
critical CTOA criterion. Rice et al. [118], extending the work of Rice and Sorensen 
[121], obtained Prandtl slip-line field solutions for a Mode I, plane-strain growing 
crack in a nonhardening elastic-plastic solid. Based on the asymptotic solutions, 
they proposed that a similar geometric profile of crack opening very near the tip 
is maintained during crack growth. The criterion for continuing crack growth in 
[121, 118] is: 


A 

d 


a dJ 
<7 () da 


+ 3 > 


eR 

~d 


( 2 . 1 ) 


where 8 C is the critical CTOD measured at a small characteristic distance d behind 
the growing crack tip 1 . The dJ/da represents the rate of external applied loading 
during crack growth. The <7o, E, and e in Equation (2.1) are the yield stress, 
the elastic modulus, and the natural logarithm base, respectively. The length 
parameter R and material parameters a and [3 are to be determined by tests or 
numerical analyses. This asymptotic field crack growth criterion is equivalent to 
the CTOA fracture criterion. 

The ability of the CTOA fracture criterion to simulate elastic-plastic crack 
growth has been verified by many numerical analyses, de Koning [36] and Anderson 
[2] were among the first to demonstrate that CTOA can be used to characterize 


x The characteristic distance is called A l and r m in [121] and [118], respectively. 
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stable crack growth. Newman [90, 91], using a 2D elastic-plastic finite element 
analysis, was able to simulate stable crack growth and predict residual strength 
within 10% for the test configurations shown in Figure 2.1. 

Recently, the CTOA fracture criterion has been extensively verified through 
tests and analyses for various loading, geometry, and crack configurations 2 . Tests 
on aluminum alloys [34, 94], 2D [93, 35, 30, 33], thin-shell [20, 23, 22], and 3D 
[27, 29, 28] elastic-plastic crack growth analyses have been conducted to assess 
the CTOA fracture criterion for aging aircraft applications. These studies will be 
discussed in Chapter 3. 

The CTOA fracture criterion is used in this study to characterize stable crack 
growth in thin-sheet metallic materials. The definition of CTOA as suggested by 
Newman [90] is adopted. For Mode-I only deformations, it is defined as (Fig- 
ure 2.2): 

CTOA = 2 tan -1 — (2.2) 

where S is the CTOD measured at a specific distance, d, behind the crack tip. For 
mixed-mode problems, the opening angle is obtained from the cross product of two 
vectors: 

CTOA = sin” 1 f** (2.3) 

IMIiIrI 

where a and b are the vectors from the crack tip to crack edges at a specific 
distance, d, behind the crack tip. 


2.4 Elastic-Plastic Crack Growth 

Stable crack growth seems to be an inherent feature of elastic-plastic materials be- 
cause of the occurrence of permanent plastic deformations during unloading [115]. 
This effect can be demonstrated by global energy dissipation or by the local resid- 
ual plastic deformations. The energy dissipation effect on stable crack growth is 
illustrated by considering the example used in [115]. Suppose two materials have 


2 Most of the test data are available from the Internet at irwin.larc.nasa.gov. 
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Figure 2.3: Illustration of crack growth in nonlinear elastic and elastic-plastic ma- 
terials. 

the same uniaxial stress-strain curves; one is an idealized nonlinear elastic mate- 
rial and the other is an elastic-plastic material. For cases without crack growth, 
the same CTOA and strain concentration will occur in the two materials as illus- 
trated in Figure 2.3, STAGE 0. As the crack propagates in the nonlinear elastic 
material, deformation fields need to be readjusted and the same crack tip open- 
ing profile would occur for the new crack tip location [115]. This is not the case 
for the elastic-plastic material because a large part of the energy is consumed by 
plastic dissipation with far less strain recovered during unloading. Thus, a smaller 
CTOA is obtained after crack growth (STAGE 1). Further increase of the applied 
loading is needed to open the crack (STAGE 2) and causes stable crack growth 
in the elastic-plastic material. Fracture instability will occur as the crack reaches 
a steady-state condition in which the crack continually advances without further 
increase in load. If the analysis is performed under displacement control, then a 
reduction in applied load is required to maintain a constant CTOA for continuous 
crack growth. Hereafter, CTOA a is the crack tip opening angle measured immedi- 
ately after propagation, STAGE 1. CTOA 6 is denoted as the increase in crack tip 
opening angle required to reach the critical value (CTOA c ). Thus, 

CTOA a + CTOA 6 = CTOA c (2.4) 

satisfies the fracture criterion for crack propagation, and the condition 

CTOA a = CTOA c (2.5) 

indicates the occurrence of fracture instability for the analysis under load control. 

Another related factor for stable crack growth is the plastic wake effect caused 
by the residual plastic deformations [90]. As the crack grows, the plastic zone 
behind the crack tip unloads to an elastic state leaving the appropriate plastic 
wake behind the advancing crack tip. This effect results in resistance to crack 
tip opening as illustrated in Figure 2.4. The dashed curves in the plastic wake 
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Figure 2.4: Illustration of plastic wake effect caused by crack growth. 

region show what the crack opening profile would have been if residual plastic 
deformations had not been retained in the material behind the advancing crack 
tip. This phenomenon is also essential for simulating the initiation of stable crack 
growth associated with high fatigue stress prior to tearing [30]. 

It is of interest to note that a different density of residual plastic deformation 
would develop as the crack propagates under plane stress compared to plane strain 
conditions. Based on the Prandtl field together with an elastic sector following the 
centered fan, Rice et al. [118, 117] have shown that asymptotic plastic strains for 
quasi-static crack growth under plane strain conditions in elastic-perfectly plastic 
materials are: 

( a/3 = F a!3{0)H^) r-> 0 (2.6) 

where (r, 0) is a local polar coordinate system with the origin at the crack tip 
and F a p(6) are functions determined from an asymptotic angular integration of 
the plastic strain rate. No complete asymptotic solutions are available for plane 
stress conditions, but only two types of plastic sectors can exist near the crack 
tip; one is the centered fan and the other is constant stress [117, 86]. For the case 
with the center fan sectors, Rice [115] shows that plastic strains under plane stress 
conditions are: 

£j = Gii{6) ln 2 (-) 0 = 0 r-+ 0 (2.7) 

where Gij are scalars from solutions with a centered fan on the 0 = 0 ray. By 
comparing Equations (2.6) and (2.7), one finds that plastic strains on the 0 = 0 
ray have a stronger singularity in plane stress than in plane strain. This observa- 
tion gives a preliminary indication that higher residual plastic deformations may 
occur under plane stress conditions leading to higher resistance to the opening of 
a growing crack. 

2.5 Link-up and Residual Strength Analysis with 
CTOA Fracture Criterion 

Since analyses based on the CTOA fracture criterion are direct simulations of real- 
istic crack growth, multiple crack growth interaction and link-up are automatically 
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Figure 2.5: Residual strength diagram under load control. 



Figure 2.6: Residual strength diagram under displacement control. 


captured as the crack propagates. Residual strength of a damaged structure is also 
obtained directly from the crack growth data. For tests conducted under load con- 
trol with the plastic zone well-confined by the elastic region, fracture instability 
is reached when no further increase of the applied load is required to maintain 
quasi-static crack extension. For tests under displacement control, the maximum 
load carrying capacity of a structure occurs followed by a reduction in load during 
continued crack growth. The residual strength diagrams corresponding to load 
control and displacement control are illustrated in Figures 2.5 and 2.6, respec- 
tively. Note that under displacement control, the load instability occurs before 
Equation (2.5) is satisfied. By comparing Figures 2.5 and 2.6, one finds that the 
load-crack extension curve up to residual strength is obtained under load control. 
On the other hand, the curve after residual strength can be obtained from the 
displacement-control technique. 
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2.6 Guidelines for Using the CTOA Fracture Cri- 
terion 

Guidelines for using the CTOA criterion to predict fracture behavior and residual 
strength of built-up aircraft fuselages are presented below 3 for future reference. 
The guidelines focus on (1) how to obtain a valid CTOA value from laboratory 
fracture tests for a given thickness of thin-sheet metallic material, and (2) how to 
correlate or fine-tune numerical analyses based on laboratory tests to predict the 
fracture behavior and residual strength of complex structures. 

1. Conduct tension tests to accurately describe the stress-strain behavior of the 
material. 

2. Conduct fracture tests on coupon specimens. Measure the CTOA c and record 
crack growth data during stable tearing. The experimental CTOA c measure- 
ments typically have a scatter band of ±0.5° to ±1.0°. The material of the 
specimen should be the same alloy, temper, and thickness as the material 
of the complex structure. The specimen should be large enough to allow 
significant crack growth prior to reaching residual strength. 

3. Use thin-shell elastic-plastic finite element analyses accounting for 3D con- 
straint effects developed at the crack tip 4 to simulate fracture behavior. Com- 
pare the predicted load versus crack growth to the experimental data and 
determine the value of CTOA c that best correlates the experimental data 
and numerical results. Note that: 

• The characteristic distance, d. used in analyses should be the same as 
the one used in experimental measurements. 

• Agreement of the CTOA c values obtained independently from exper- 
imental measurements and numerical analyses would greatly increase 
the confidence in the chosen CTOA c value. 

4. Create a finite element model for the structure to be analyzed. Use the 
previously determined CTOA c value to predict the fracture behavior. The 
model should have the same size and type of crack tip elements as the one 
used in the coupon test correlation. 

The validity of these guidelines as applied to fuselage structures will be examined 
extensively through Chapters 3 and 4. 

3 These guidelines follow closely the recommendations made in [33, 29]. 

4 One way to consider the 3D constraint effects in a thin shell analysis is to use 
the plane strain core concept (see Figure 3.2). 
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2.7 Summary 

Various LEFM and EPFM versions of fracture mechanics methods to characterize 
fracturing processes in thin-sheet metals were reviewed and critiqued in Section 2.2. 
Among them, the CTOA fracture criterion with a 3D elastic-plastic analysis was 
found to be a superior one because of its relative independence of the geometry of 
the structure, the length of the crack, and the presence of multiple cracks. 

Previous experimental, analytical, and numerical studies for the CTOA fracture 
criterion were reviewed in Section 2.3. The definition of the CTOA was given in 
Equation (2.2) for Mode-I only deformations and in Equation (2.3) for general 
mixed-mode problems. 

The CTOA-driven elastic-plastic crack growth was studied in Section 2.4. The 
inherent feature of stable crack growth in elastic-plastic materials was discussed 
using the energy dissipation and residual plastic deformation. Residual plastic 
strains from elastic-plastic crack growth under plane stress and plane strain con- 
ditions were studied using the Prandtl field. Higher residual plastic deformations 
were found under plane stress conditions. As a result, a higher fracture resistance 
of a growing crack may occur under plane stress conditions. This behavior will be 
further examined numerically in Chapter 3. 

Analyses of link-up of multiple cracks and residual strength of damaged struc- 
tures using the CTOA fracture criterion were discussed in Section 2.5. Finally, 
guidelines for using the CTOA criterion calibrated from coupon tests to predict 
fracture behavior of built-up aircraft fuselages were presented in Section 2.6. 



Chapter 3 

Residual Strength Analysis of a 
Flat Panel with Self-Similar 
Elastic-Plastic Crack Growth 


Elastic-plastic crack growth simulations and residual strength prediction of flat 
panel coupon tests are studied in this chapter. The purposes of this chapter are 
to: 

1. further review and discuss some recent activities of using the crack tip open- 
ing angle (CTOA) fracture criterion for aging aircraft applications, 

2. model the fracturing processes in middle-crack tension (MT) specimens using 
elastic-plastic, thin shell finite element analyses, 

3. explore the need to incorporate the three-dimensional constraint effect to 
characterize fracture behavior of thin-sheet metals, and 

4. model the fracturing processes in thin-sheet specimens with multi-site dam- 
age (MSD). 

3.1 Introduction 

Tests and numerical simulations have been performed to assess the CTOA fracture 
criterion for predicting residual strength of aging aircraft. Laboratory tests were 
conducted on flat panels made of aluminum alloys [34, 94], Numerical simulations 
were conducted using two-dimensional (2D) [93, 35, 94, 30, 33], thin-shell [20, 
23, 22], and three-dimensional (3D) [27, 29, 28] finite element elastic-plastic crack 
growth analyses. We review these activities in a somewhat chronological order and 
highlight the important findings of these studies below. The latest results are used 
as a starting point for subsequent simulations in this study. 

A series of fracture tests have been conducted using a 2024-T3 aluminum al- 
loy for MT, CT, blunt notch, THCT and MSD specimens. Newman et al. [93] 
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conducted tests on 0.05, 0.07, and 0.09 inch thick, 3.0 and 11.8 inch wide MT 
and blunt notch specimens as well as 0.09 inch thick, 10 inch wide THCT spec- 
imens. The blunt-notch specimen is similar to the MT specimen except that a 
small hole is drilled at both ends of the saw cut. It is intended to assess the suit- 
ability of elastic-plastic finite element analyses with the small-strain assumption 
to model large-scale plasticity deformations. A good agreement between predicted 
and measured load versus notch-tip displacements substantiates the assumption. 

The critical values of CTOA (CTOA c ) were measured for the MT and THCT 
specimens to assert the specimen configuration independence of the fracture cri- 
terion. The THCT specimen has a stress intensity factor solution like that of a 
cracked, stiffened panel [91]. The measured CTOA c values showed higher angles at 
crack initiation, but reached the same constant value after a small transition pe- 
riod of crack growth. The agreement of CTOA c between MT and THCT specimens 
indicates that the CTOA fracture criterion is independent of specimen configura- 
tion; this was further confirmed by a follow-up study with measurements from CT 
specimens [35]. 

A 2D elastic-plastic finite element code, ZIP2D [88], and a 6.1 degree CTOA c , 
computed at 0.01875 inch behind the crack tip, were used to simulate fracture 
behavior of the MT specimens [93]. To model fatigue pre-cracking, cyclic loading- 
simulation was conducted prior to stable tearing analyses. Experimental and pre- 
dicted results showed that a higher applied stress during the fatigue tests increased 
the resistance of stable crack growth initiation. Predicted residual strengths under 
plane stress conditions were within 4% of experimental results for 3.0 and 11.8 inch 
wide MT specimens. Yet the plane stress analyses over-predicted crack extensions 
prior to limit load while the plane strain analyses under-predicted crack extensions. 

The above studies raised two important questions: 

1. What is the governing mechanism that causes higher CTOA c values during 
crack initiation? 

2. What is the governing mechanism that causes the discrepancy between 2D 
predictions and test results? 

Dawicke and Sutton [30] examined the higher values of measured CTOA c ob- 
served during crack initiation, i.e., question 1. Two independent techniques, op- 
tical microscopy (OM) and digital image correlation (DIC) were used to measure 
surface CTOA c during crack growth. The results of the two methods agreed very 
well. Fatigue marker loads and a scanning electron microscope were used to exam- 
ine the fracture morphology and sequences of crack front profiles. For specimens 
under low magnitude of fatigue stress prior to tearing, crack surfaces underwent a 
transition from flat-to-slant crack growth. A schematic of the transition is shown 
in Figure 3.1. During the transition period, the CTOA c values were high and sig- 
nificant tunneling occurred. After an amount of crack growth equal to about the 
specimen thickness, CTOA c reached a constant value. After crack growth equal 
to about twice the thickness, crack tunneling stabilized. For specimens that were 
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Figure 3.1: Schematic of fracture surface indicating transition from a flat to a slant 
crack plane (after [93]). 

pre-cracked under a high magnitude of fatigue stress, a 45-degree, slant, through- 
thickness initial crack was formed prior to tearing. During the crack initiation 
period, the CTOA c values of specimens with high fatigue stress were lower than 
the ones with low fatigue stress. But the same constant CTOA c value was observed 
after crack growth equal to about the specimen thickness. 

The discrepancy between 2D predictions and test results, i.e., question 2, was 
thought to be related to the 3D constraint effect. Although thin-sheet structures 
behave essentially in plane stress, the constraint due to the finite thickness of the 
specimens can cause the regions local to the crack tip to approach plane strain 
conditions [56]. 

To investigate the constraint effect, 2D and 3D analyses were conducted. In 
the 2D analyses, a core of elements above and below the crack path were assigned 
as plane strain while all other elements were assigned as plane stress. The plane 
strain core concept is illustrated in Figure 3.2. 

In their early attempt, Dawicke et al. [35, 94] used 2D finite element analyses 
with a 6.0 degree CTOA c computed at 0.02 inch behind the crack tip and a plane 
strain core height equal to 0.2 inch to simulate fracture behavior with the constraint 
effect. They showed that the use of a plane strain core was essential to accurately 
model crack growth. The predicted residual strengths were within 2% for 3 and 12 
inch wide, 0.09 inch thick MT specimens and within 4% for 6 inch wide, 0.09 inch 
CT specimens. For 20 inch wide, 0.04 inch thick MSD specimens, 2D analyses with 
a 5.1 degree CTOA c showed excellent agreement of link-up and residual strength 
between predictions [94] and test results [13]. 

Dawicke et al. [27, 29] further studied the constraint effect using 3D finite 
element analyses with a 5.25 degree CTOA c computed at 0.04 inch behind the 
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crack tip. The 3D analyses successfully simulated fracture behavior of 2.0, 4.0, 
6.0, and 8.0 inch wide CT specimens, 1.2, 3.0, 6.0, 12.0, 24.0, and 60.0 inch wide 
MT specimens, and 12.0 inch wide MSD specimens made of 0.09 inch thick, 2024- 
T3 aluminum alloy. A plane strain core height of 0.12 inch was required for 2D 
analyses to match the measured results and the 3D fracture predictions. 

In the following, the MT and MSD tests are studied. The FRANC3D/STAGS 
program is used to simulate fracture behavior and predict residual strength using 
the guidelines derived from the above 2D and 3D studies. 

3.2 Experimental Procedures and Test Configu- 
rations 

Fracture tests of MT specimens were conducted by the Mechanics of Materials 
Branch at NASA Langley Research Center [34, 27, 29]. The test specimens were 
made of 0.09 inch thick 2024-T3 aluminum alloy. All specimens were fatigue pre- 
cracked in the L-T orientation with a low stress level that results in a stress inten- 
sity factor range of AK = 7 ksiVinch. For specimens with a single crack, different 
widths of panels equal to 3 inch, 12 inch, and 24 inch with a crack-length to width 
ratio equal to 1/3 were tested (Figure 3.3). For cases with multiple cracks, only 
the 12 inch wide specimens with two to five near collinear cracks as illustrated in 
Figure 3.4 were tested. All tests were conducted under displacement control with 
guide plates to prevent out-of-plane buckling. Both OM and DIC techniques were 
used to measure the CTOA c during stable crack growth [34] . Results for MT and 
CT specimens are shown in Figure 3.5. The CTOA c rapidly reaches a constant 
value with a scatter band about ±1.0°. 
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Figure 3.3: Test configurations of MT specimens. 
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Figure 3.4: Test configurations of 12 inch wide specimens with multiple cracks. 
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Figure 3.5: Surface measurements of CTOA c (after [29]). 

3.3 Numerical Simulations of MT Specimens 

Fracture processes of MT specimens were simulated first. To investigate panel 
size effects, numerical simulations of 60 inch wide panels with the same crack- 
length-to-width ratio were also performed. Elastic-plastic finite element analyses 
based on incremental flow theory with the von Mises yield criterion and the small 
strain assumption were used to capture the active plastic zone and the plastic wake 
during stable crack propagation. A piecewise linear representation was used for 
the uniaxial stress-strain curve for 2024-T3 aluminum (Figure 3.6). The CTOA c 
used in this study was 5.25 degrees measured 0.04 inch behind the crack tip. This 
particular CTOA value was provided by Dawicke and Newman [27, 29] based on 
3D simulations of CT specimens 1 . Upon satisfaction of the fracture criterion, 
nodal release and load (or displacement) relaxation techniques were employed to 
propagate the crack. Because of the double symmetry of the geometry and loading, 
only one quarter of the specimen with imposed symmetry boundary conditions 
was modeled. Out-of-plane displacements were suppressed. Displacement-based 
four-noded and five-noded quadrilateral shell elements having C 1 continuity were 
used [112]. These elements are intended to model thin shell structures for which 
transverse shear deformation is not important. Each node of the element has six 
degrees of freedom including three translations and three rotations. A special five- 
noded shell element, formulated by combining two four-noded elements and using 
linear constraint along the edge to eliminate the dependent node, was used to 

x As noted by Dawicke and Newman [29], the fracture behavior of the CT spec- 
imen is more sensitive to small changes in CTOA c than the MT specimen; thus 
the CT specimen is more suitable to correlate predicted and measured results. 
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Figure 3.6: Piecewise linear representation of the uniaxial stress-strain curve for 
2024-T3 aluminum. 

transition from locally refined zones around the crack path to a coarse mesh away 
from the crack. 

A convergence study was conducted to determine the sensitivity of the predicted 
residual strength to the element size along the crack extension path. Three meshes 
for the 24 inch wide panel were created with crack tip element sizes of 0.04 inch, 
0.02 inch, and 0.01 inch. For all crack growth and residual strength analyses, 
the CTOA was evaluated at 0.04 inch behind the crack tip to be consistent with 
experimental measurements. A finite element mesh with 0.04 inch square crack 
tip elements for the 24 inch wide panel is shown in Figure 3.7. Predicted crack 
growth results for cases with 0.04 inch and 0.02 inch crack tip elements as well as 
predicted residual strengths for all three cases are shown in Figure 3.8. Although 
some discrepancy was observed at the early stage of stable tearing, the predicted 
results exhibited little influence of mesh size after a relatively small amount of 
stable crack growth. More importantly, the predicted residual strength was very 
insensitive to crack tip element size. Thus, all the remaining meshes used in this 
study had 0.04 inch crack tip elements. 

3.3.1 Numerical Results 

Figure 3.9 shows two predicted crack opening profiles for the 24 inch wide panel. 
The angles were computed immediately after propagation (i.e., CTOA a , see Fig- 
ure 2.3) with relaxation procedures completed and before increasing the applied 
displacement. The two CTOA a values correspond to (1) the angle after the first 
increment of crack growth, and (2) the angle after the specimen reaches its residual 
strength. As shown in the figure, CTOA„ is much smaller than the critical angle 
after the first crack growth increment. This clearly demonstrates the permanent 
plastic deformation effects on stable crack growth in the elastic-plastic material 
(c/. section 2.4). As the crack propagates, CTOA a increases. Since the analyses 
were conducted under displacement control, the CTOA a at residual strength is less 
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Figure 3.7: Finite element mesh for 24 inch wide MT specimen and detail along 
crack path. 




Figure 3.8: Convergence study: predicted crack growth and predicted residual 
strength for 24 inch wide panel with different crack tip element sizes. 
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Figure 3.9: Crack opening profiles and CTOA a after the first crack growth incre- 
ment and after reaching the residual strength for 24 inch wide panels. 

than, but approaching its critical value. 

Comparisons between numerical results and experimental measurements for 
the applied stress versus half crack extension are shown in Figure 3.10. Results of 
predicted residual strength are comparable to experimental measurements, but as 
the width of the panel increases, the relative difference between experimental mea- 
surements and numerical predictions increases. Figure 3.11 depicts the predicted 
plastic zone as the specimens reach their ultimate strength. Two distinct phenom- 
ena are observed. For small specimens, plastic zones reach the free edge and the 
limit load is attained due to net section yielding. In contrast, for large specimens, 
plastic zones are well-confined by the elastic region and residual strength is reached 
near the fracture instability of the specimens. 

3.3.2 Discussion 

As shown in Figure 3.10, the relative difference in residual strength between exper- 
imental and numerical results increases as the width of the panel increases. This 
discrepancy is believed to be due to the three-dimensional nature of the stresses 
around the crack tip, a result of constraint effects due to the finite thickness of 
the panels [56, 31]. Numerical results using plane strain, plane stress with a plane 
strain core height (see Figure 3.2) equal to 0.12 inch, and three-dimensional finite 
element analyses obtained from [27, 29] were studied to further demonstrate con- 
straint effects on residual strength predictions. Predicted results shown in Table 3.1 
and Figure 3.12 suggest that: 

• thin shell finite element analysis, behaving essentially in plane stress, tends 
to over-predict the residual strength as the width of the panel increases; 

• plane strain analysis over-predicts the residual strength of small specimens, 
but under-estimates it for large specimens; 



applied stress (ksi) applied stress (ksi) 
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half crack extension, Aa (in.) half crack extension, Aa (in.) 

Figure 3.10: Comparisons between experimental measurements and numerical pre- 
dictions of applied stress versus half crack extension for various sizes 
of specimens. 
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Table 3.1: Comparisons of Residual Strength Predictions of MT Specimens 
(unit: ksi) 
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Figure 3.12: Predicted results of thin shell, plane strain, plane stress with a plane 
strain core, and 3D analyses compared with experimental measure- 
ments. 

• 2D plane stress analysis with a plane strain core and 3D analysis properly 
account for constraint effects. The predicted results follow the trend of ex- 
perimental measurements even for wide panels. 

The cross over between plane stress and plane strain in predicting residual 
strength as the specimen size increases is an interesting topic. Based on the pre- 
dicted plasticity distribution in Figure 3.11, the net section yielding mechanism 
seems to dominate the residual strength prediction of small specimens. This may 
explain why the plane strain analysis predicts a higher residual strength for small 
specimens because the effective yield stress in plane strain is larger than that in 
plane stress. Thus, a further increase of remote stresses under plane strain con- 
ditions is needed for specimens to reach the point of net section yielding. For 
larger specimens, residual strength is governed by stable crack growth and frac- 
ture. As one would expect from the thickness effects on K c in LEFM [9], materials 
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in plane stress have higher fracture toughness than materials in plane strain. Re- 
cent micromechanics-based, 3D analysis of ductile crack growth in a thin plate 
with a Gurson-type model also showed that, although the crack growth resistance 
at first increases with increasing plate thickness, the resistance to crack growth 
decreases after a small amount of crack extension [84]. 

For CTOA-driven ductile crack growth, stresses and strains under plane stress 
and plane strain conditions have not been studied in sufficient detail to clarify the 
issue. A possible cause of higher crack growth resistance in plane stress may be 
related to the residual plastic deformation effects. Based on asymptotic solutions 
for cracks growing in an incompressible elastic-perfectly plastic material under 
Mode I loading (Equations (2.6) and (2.7)), larger residual plastic deformations 
would occur under plane stress than plane strain conditions leading to higher crack 
growth resistance. 


3.4 Numerical Simulations of Specimens with Mul- 
tiple Cracks 

Numerical simulations of tests with multiple cracks using the CTOA fracture cri- 
terion are straightforward extensions of single crack specimen simulations. The 
same fracture criterion (CTOA c = 5.25 degrees measured 0.04 inch behind the 
crack tip) was used to simulate stable crack growth and the link-up of multiple 
cracks, and to predict the residual strength. No supplementary criterion is needed. 
Multiple crack test configurations as shown in Figure 3.4 were modeled and the 
fracture processes were simulated. Note that the symmetry conditions along the 
vertical central line of the specimens (see Figure 3.4) are no longer valid due to 
the various lengths of fatigue pre-cracks; thus, at least one half of the specimen 
needs to be modeled. A finite element mesh for test configuration b is shown in 
Figure 3.13. Mesh patterns around the multiple cracks are similar to those of the 
single crack models. 

3.4.1 Numerical Results and Discussion 

Numerical results and experimental measurements for the applied stress versus 
half crack extension for test configuration b and d are shown in Figure 3.14. Two 
distinct applied load versus crack growth history curves are predicted. For test 
configuration a, b, and c, link-up of cracks happens before the specimens reach their 
residual strength. For test configurations d and e, the limit load is attained before 
link-up. These numerical predictions agree with observations from the fracture 
tests. 

Again, plastic deformation plays an important role in the fracture process. 
Figure 3.15 shows the plastic zone evolution of test configuration b during stable 
crack growth. The inherent residual plastic deformations during crack growth are 
clearly demonstrated through the deformed shapes. Figure 3.16 summarizes the 
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Figure 3.13: Finite element mesh for the test configuration b (12 inch wide speci- 
men with two cracks). 



crack extension, Aa (in. ) crack extension, Aa (in.) 


Figure 3.14: Predicted applied stress versus crack extension for test configuration 
b and d. 







35 



(1) first increment (2) before link-up 



(3) after link-up 



Figure 3.15: Crack opening profile(s) and plastic zone evolution of test configu- 
ration b during crack growth: (1) at the first increment, (2) before 
link-up, (3) after link-up, and (4) reaching the residual strength. 


relative difference between predicted results and experimental measurements. The 
predicted residual strength of all five MSD simulations agrees very well (within 3%) 
with experimental data. The predicted link-up load is comparable to experimental 
measurements, but the difference is larger than that for the residual strength. 
Reasons for the discrepancy may be related to the difficulty in measuring link-up 
load during the fracture tests. 

It is of practical importance to characterize the reduction in residual strength 
caused by MSD [50]. Figure 3.17 plots numerical predictions of residual strength 
versus lead crack length for cases with and without small cracks. A loss of residual 
strength due to the presence of multiple small cracks is observed. 


3.5 Summary 

Stable crack growth and residual strength prediction for the flat panel tests were 
performed. The CTOA criterion was used to characterize the elastic-plastic crack 
growth in thin-sheet metals. The major findings of the flat panel study are: 

1. Two distinct failure mechanisms are observed for MT specimens. For small 
specimens, plastic zones reach the free boundary and the limit load is attained 
due to net section yielding. For large specimens, plastic zones are well- 
confined by the elastic region and residual strength is reached due to the 
fracture instability of the specimens. 

2. Constraint effects caused by the finite thickness of the plates provide a reason- 
able explanation for the increase of the relative difference between predicted 
residual strength from thin shell analyses and experimental measurements as 
the panel size increases. 



Residual Strength Difference (%) 


36 



Figure 3.16: Relative difference of residual strength and link-up load between pre- 
dicted results and experimental measurements for specimens with 
multiple cracks. 
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Figure 3.17: Loss of residual strength due to the presence of small cracks. 
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3. Predicted link-up load and residual strength are in good agreement with 
experimental measurements for the panels with multiple cracks. A loss of 
residual strength due to the presence of multiple small cracks is observed. 

4. The CTOA fracture criterion, combined with elastic-plastic shell finite ele- 
ment analyses, proves to be a rational and rigorous simulation tool to charac- 
terize stable crack growth and to predict the residual strength of flat panels 
with single and multiple cracks. 



Chapter 4 

Residual Strength Analysis of 
Fuselage Structures with 
Self-Similar Crack Growth 


In this chapter, the crack tip opening angle (CTOA) fracture criterion obtained 
from coupon tests is used to predict fracture behavior and residual strength of 
built-up aircraft fuselages that are subjected to widespread fatigue damage (WFD). 
Two fuselage models are investigated. The first example is a generic narrow body, 
lap-jointed fuselage with stringers and frames but without tear straps. This rel- 
atively simple, built-up configuration is used to demonstrate the ability of the 
FRANC3D/STAGS system to predict residual strength of fuselage structures sub- 
jected to WFD. The second example is a detailed analysis of a wide body, lap- 
jointed fuselage panel with tear straps, stringers, stringer clips, and frames. The 
analyses focus on simulations of single crack growth and multi-site damage (MSD) 
in a fuselage panel conducted in a full-scale, wide body, pressure test fixture [49, 50]. 
This example is intended to validate the analysis methodology by directly com- 
paring numerical predictions with experimental measurements on actual fuselage 
structures. 


4.1 Demonstration Example: A Generic Narrow 
Body Fuselage Panel 

A relatively simple built-up narrow body fuselage configuration was modeled. The 
example demonstrates the analysis capability to predict the residual strength of a 
pressurized fuselage, subjected to WFD and corrosion damage [25, 24], The prob- 
lem chosen for analysis was a three stringer wide, three frame long fuselage panel. 
The panel section had a radius of curvature of 72 inches. It contained a lap joint at 
the central stringer. The lap joint was a typical three row configuration with 3/16 
inch diameter countersunk-head rivets. The other two stringers were spot-welded 
to the skin. The upper and lower skins were made of 0.04 inch thick, 2024-T3 
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Figure 4.1: Dimensions of generic narrow body panel for demonstration example. 

aluminum alloy. The stringers and frames were made of 7075-T6 aluminum alloy. 
Frames were simply connected to stringers by rivets. The panel configurations are 
shown in Figures 4.1 and 4.2. The frame and stringer dimensions are shown in 
Figure 4.3. 

4.1.1 Numerical Model 

All structural components including skins, stringers, and frames were modeled by 
shell elements. Each node of a shell element has six degrees of freedom. A piecewise 
linear representation was used for the uniaxial stress-strain curves for 2024-T3 
and 7075-T6 aluminum alloys (see Figures 4.4 and 4.5). Symmetric boundary 
conditions were imposed on all the boundary edges to simulate a cylinder-like 
fuselage structure. Pressure loading was applied on all the external skins. 

Both geometric and material nonlinearities were included in the analysis. The 
former captures the out-of-plane bulging deformation and the latter captures the 
active plastic zone and the plastic wake during stable crack propagation. The 
nonlinear solution algorithm consists of Newton’s method. Large rotations were 
included in the nonlinear solution by a co-rotation algorithm applied at the element 
level [96]. The Riks arc-length path following method was used to trace a solution 
past the limit points of a nonlinear response [122, 113]. 

Rivets were modeled by elastic-plastic spring elements that connect finite ele- 
ment nodes in the upper and lower skins. Each rivet was modeled with six degrees 
of freedom, corresponding to extension, shearing, bending and twisting of the rivet. 
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Figure 4.2: Detailed rivet spacing for demonstration example. 



Figure 4.3: Dimensions of stringer and frame for demonstration example (dimen- 
sions in inches, modified after [25]). 
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Figure 4.4: Demonstration example: piecewise linear representation of the uniaxial 
stress-strain curve for 2024-T3 aluminum. 
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Figure 4.5: Demonstration example: piecewise linear representation of the uniaxial 
stress-strain curve for 7075-T6 aluminum. 
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Figure 4.6: Demonstration example: rivet shear stiffness and strength. 


The stiffness of each degree of freedom was defined by prescribing a force-deflection 
curve. The axial, flexural, and torsional stiffnesses of the spring element were com- 
puted by assuming that the rivet behaves like a simple elastic rod with a diameter 
of 3/16 inch. The elastic shear stiffness of the rivet was computed by the following 
empirical relation developed by Swift [136]: 


K 


rivet 


ED 

[A + ( (D/t\ + D/l> 


(4.1) 


where E is the elastic modulus of the sheet material, D is the rivet diameter, t\ and 
t 2 are the thicknesses of the joined sheets, and A = 5.0 and C = 0.8 for aluminum 
rivets. The initial shear yielding and ultimate shear strength of the rivets were 
assumed to occur at load levels of 510 lb and 725 lb, respectively. Once a rivet 
reaches its ultimate strength, it will break and lose its load carrying capacity. The 
force-deflection curve shown in Figure 4.6 for shearing is intended to represent 
empirically the net shear stiffness of a rivet-joined sheet connection, accounting 
for bearing deformations and local yielding around the rivet [136, 150]. 

The critical crack tip opening angle (CTOA c ) was used to characterize elastic- 
plastic crack growth and to predict residual strength. The CTOA c used in this 
example was 5.7 degrees measured 0.04 inch behind the crack tip with a plane 
strain core height equal to 0.08 inch 1 . Six different crack configurations with 
various lengths of lead and MSD cracks were studied. The initial configurations 
prior to crack growth were: 


1. a 7.14-inch lead crack, 


2. a 7.14-inch lead crack with 0.025 inch MSD cracks emanating from both sides 
of a fastener hole, 

1 Since no experimental crack growth data were available, this particular CTOA c 
value was estimated based on the 5.25 degrees used in 0.09 inch thick, 2024-T3 
bare material in Chapter 3. The plane strain core height was assumed to be twice 
the sheet thickness. 
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Figure 4.7: Crack configurations with a 10-inch initial lead crack and MSD (exter- 
nal view). 

3. a 7.14-inch lead crack with 0.046 inch MSD cracks emanating from both sides 
of a fastener hole, 

4. a 10-inch lead crack, 

5. a 10-inch lead crack with 0.025 inch MSD cracks emanating from both sides 
of a fastener hole, and 

6. a 10-inch lead crack with 0.046 inch MSD cracks emanating from both sides 
of a fastener hole. 

The lead crack was located symmetrically about the central frame line. The MSD 
pattern was symmetric about the lead crack at the 3 rivets in front of the lead 
crack. The lead and MSD cracks were located along the upper rivet row in the 
upper skin of the joint. The crack configurations with a 10-inch initial lead crack 
are shown in Figure 4.7. Since rivet holes were not modeled explicitly in the finite 
element model, a small crack with a length equal to the rivet diameter plus the 
MSD length was used to model the MSD crack. 

A mesh pattern with 0.04 inch crack tip elements was used. This pattern is 
similar to the one used in the flat panel simulation (c/. Figure 3.7). A finite element 
mesh for the model is shown in Figures 4.8 and 4.9. In addition to the effects of 
WFD, material thinning due to corrosion damage was also studied. The effect of 
material thinning was modeled by a uniform reduction in thickness of the upper 
skin at the lap joint in the two center bays. 
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Figure 4.8: Finite element mesh for demonstration example. 


initial crack prescribed crack growth path 




Figure 4.9: Detailed mesh around crack path for demonstration example. 
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4.1.2 Numerical Results and Discussion 

Figure 4.10 shows the predicted results of the operating pressure loading versus the 
total crack extension for all the cases conducted in this study. Predicted residual 
strengths summarized in Figure 4.11 indicate: 

• The MSD cracks significantly reduce the residual strength of the fuselage 
panel. A 21.8 to 28.0% loss of residual strength due to the presence of small 
MSD is observed. 

• A 10% uniform thickness degradation due to corrosion damage reduces the 
residual strength by 3.4 to 9.0%. The coupling of WFD and corrosion damage 
leads to the most severe damage scenario. 

• In general, increasing the lead and MSD crack lengths reduces the residual 
strengths. However, for the cases with a 10-inch initial lead crack, residual 
strength seems to be relatively insensitive to the MSD crack sizes. 

The deformed structure at residual strength for the case with a 10-inch initial 
lead crack but without MSD and corrosion damage is shown in Figure 4.12. Out- 
of-plane bulging is observed in the skin crack edges. Because of the stiffness of 
the stringer, the bulging at the lower crack edge is much smaller than the oppos- 
ing edge. The unsymmetric out-of-plane bulging thus leads to an anti-symmetric 
bending deformation field at the crack tips [108]. 

Figures 4.13 and 4.14 depict the predicted plastic zones for the cases with a 
10-inch initial lead crack as the panel reaches its residual strength. As shown 
in Figure 4.13, the evolving plastic zones are well-confined by the elastic regions 
within the frames. For the case without MSD, dominant plastic zones accom- 
panying the lead crack tips are observed. For the case with MSD, plastic zones 
are developed at the multiple crack tips. The plasticity distributions are highly 
influenced by the multiple crack interactions. 


4.2 Validation Example: A Generic Wide Body 
Fuselage Panel 

Full-scale pressurized panel tests described in [49, 50] were simulated. The tests, 
funded by the Federal Aviation Administration (FA A), and performed by the Boe- 
ing Commercial Airplane Group, were intended to characterize crack growth in a 
generic wide body, lap-jointed fuselage configuration, subjected to WFD. Detailed 
analyses using the FRANC3D/STAGS program were conducted to validate the 
analysis methodology. Stress distributions were compared with strain gage read- 
ings. Predicted stable crack growth and residual strength results were compared 
with experimental measurements. 



Operating Pressure (psi) Operating Pressure (psi) 
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— □ — No MSD 

— A — 0.025 in. MSD 

— O — 0.046 in. MSD 



Total Crack Extension (inch) Total Crack Extension (inch) 


(a) (b) 



Total Crack Extension (inch) Total Crack Extension (inch) 


(c) (d) 

Figure 4.10: Predicted operating pressure versus total crack extension for the 
demonstration example: (a) 7.14-inch initial lead crack, (b) 7.14-inch 
initial lead crack with corrosion damage, (c) 10-inch initial lead crack, 
and (d) 10-inch initial lead crack with corrosion damage. 
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Initial Lead Crack Length (inch) 


□ No MSD 
A 0.025 in. MSD 
O 0.046 in. MSD 
a No MSD; corrosion 
^ 0.025 in. MSD; corrosion 

• 0.046 in. MSD; corrosion 

curve fit 


Figure 4.11: Predicted residual strength versus initial lead crack length. 



Figure 4.12: Deformed shape of the demonstration example (pressure = 15.3 psi, 
magnification factor = 5.0). 
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Effective Stress (ksi) 




60.0 (C u ) 

53.0 

46.0 

39.0 (o y ) 


Figure 4.13: Predicted plastic zones for 10-inch initial lead crack without MSD 
(pressure = 15.3 psi, magnification factor = 5.0). 


Effective Stress (ksi) 

60.0 (O u ) 

53.0 

46.0 

39.0 (O y ) 

Figure 4.14: Predicted plastic zones for 10-inch initial lead crack with 0.025 inch 
MSD (pressure = 11.3 psi, magnification factor = 5.0). 
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Figure 4.15: Generic narrow and wide body test fixtures (after [82]). 

4.2.1 Full-Scale Fuselage Panel Testing 

Two generic pressure test fixtures were fabricated by the Boeing Commercial Air- 
plane Group. One fixture had a radius of curvature of 74 inches to match narrow 
body airplanes and the other had a radius of curvature of 127 inches to match wide 
body airplanes. The test fixtures are shown in Figure 4.15. One end of each fixture 
was mounted in a rigid framework and the other on rollers to allow longitudinal 
displacement. Removable test sections were inserted in cutouts in the fixtures. 
Tests were conducted under pressure loading only, using air as the pressurizing 
medium. The full-scale fuselage panel tests investigated in this section were per- 
formed on the wide body pressure test fixture. A brief overview of the panel tests 
is described below. More information about the fixtures and tests can be found in 
[82, 85, 49, 50], 

Two identical curved lap-jointed panels were fabricated. The test panels were 
designed to simulate typical wide body fuselage crown structures consisting of 
bonded tear straps and floating frames connected to hat section stringers with 
stringer clips. Skins and tear straps were made of 0.063 inch thick, 2024-T3 clad 
aluminum alloy. Stringers, frames, and stringer clips were made of 7075-T6 clad 
aluminum alloy. The skins were joined by the lap joints. The joint was a typ- 
ical three row configuration assembled using standard 3/16 inch diameter, 100° 
countersunk-head rivets. The tear straps were hot-bonded to the skins at each 
frame station. The outer and inner tear straps were overlapped above the lap 
joint. The dimensions of the panels are shown in Figures 4.16, 4.17, and 4.18. The 
dimensions of frames, stringers, and stringer clips are shown in Figures 4.19 and 
4.20. 

A five-inch initial saw cut was inserted along the upper rivet row in the outer 
skin. For the panel with MSD cracks, small sawcuts were inserted in the outer 
skin after the rivet holes had been drilled, but prior to the application of the fay 
sealant and rivet installation. The panels were subjected to pressure cycling until 
the length of the crack reached about two frame bays. The central frame was 
then cut and the residual strength tests were conducted. Rosette strain gages were 
installed back-to-back on the skins and tear straps in the vicinity of the lap joint. 
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Figure 4.16: Validation example: test panel dimensions. 




Figure 4.17: Validation example: detailed panel dimensions (after [49]). 
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Fay Surface Sealant 
(Typ both lap joints) 



1.31 inches 


1.35 inches J 


1.4- 1.5d Minimum 
button dia. on lap joint 



v I 1YJA16N hj 
t 0.4 0 inch k . c ' \ 


2 Eq. 
Spaces 


R 0.38\ 
inch 
(Typ) 


0.44 

inch 


4.20 inches 
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2.60 ' 


Skin 




YJA 

6N 

C 



inches 


XF 

6N 

C 



BACB30IMW6K (6AI-4V Hi-Lok 1 00° CSK Head) 3/1 6 inch diameter 
install per 69B15265 method “A" 

BACR15CE6D (2017— T4 Rivet 100° CSK Head) 3/16 inch diameter 


Figure 4.18: Validation example: detailed rivet spacing (after [49]). 
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Frame line 



Figure 4.19: Validation example: dimensions of frame and stringer clip (dimensions 
in inches). 



h — 1.0 — H 


Figure 4.20: Validation example: dimensions of stringer (dimensions in inches). 
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4.2.2 Numerical Model 

All structural components including skins, stringers, and frames were modeled 
by displacement-based four-noded or five-noded quadrilateral shell elements [113, 
112]. To analyze the panel tests with reasonable computer resources and sufficient 
accuracy, a global-local approach was used. Figure 4.21 shows the typical finite el- 
ement meshes for the two hierarchical modeling levels employed in the simulations. 
A 12-stringer-bay wide and 5-frame-bay long panel, which is about the size of the 
test panel, was modeled at the global level. A lxl bay stiffened panel was modeled 
at the local level. The local model differed from the global model in the finite 
element mesh density and the detailed geometric modeling of the cross sectional 
shapes of stringers and frames. 

Pressure loading was applied on all the external skins. Symmetric boundary 
conditions were imposed on all the boundary edges of the global model to simulate 
a cylinder-like fuselage structure. Uniform axial expansion was allowed at one 
longitudinal end. On this boundary edge, an axial force equal to (PR/ 2) ■ L was 
assigned where P is the applied pressure, R is the radius of the panel, and L is 
the arc-length of the edge. The kinematic boundary conditions (displacements and 
rotations) applied along the boundaries of the local model were extracted from the 
global model results. In addition to these kinematic constraints, the local model 
was also subjected to internal pressure. 

A piecewise linear representation was used for the uniaxial stress-strain curves 
for 2024-T3 and 7075-T6 aluminum alloys (see Figures 4.22 and 4.23). Similar to 
the demonstration example, rivets were modeled by elastic-plastic spring elements. 
The shear force-deflection curve for the rivet is shown in Figure 4.24. Since no 
special adhesive elements were available in the STAGS element library, the adhesive 
bond between skin and tear strap was also modeled with spring elements. The shear 
stiffness for the springs was computed based on an effective area of the adhesive 
with [128]: 


K 


adhesive — 


A e ff 

ta/G a + (3/8) (t\/G + p /G)) 


(4.2) 


where A e ff is the bond area being lumped at the finite element nodal connection, 
G is the elastic shear modulus of the sheet material, G a is the elastic shear modulus 
of the adhesive, t\ and f 2 are the thicknesses of the bonded sheets, and t a is the 
thickness of the adhesive bond. Because no adhesive tests were conducted, the 
material properties of adhesive, G a and t a , were obtained from the experimental 
results in [135]. The maximum shear deflection of the adhesive bond was assumed 
to be 0.001 inch. Similar to the rivet spring, once the adhesive spring reaches 
its ultimate strength, it will break and lose its load carrying capacity. The force- 
deflection curve for shearing is shown in Figure 4.25. The axial stiffness of the 
adhesive spring was derived from the shear stiffness. The torsional and flexural 
stiffnesses of adhesive were assumed to be negligible. 

Both geometric and material nonlinearities were used in the analysis at the 
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Figure 4.22: Validation example: piecewise linear representation of the uniaxial 
stress-strain curve for 2024-T3 aluminum. 
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Figure 4.23: Validation example: piecewise linear representation of the uniaxial 
stress-strain curve for 7075-T6 aluminum. 



0 0.01 0.02 0.03 0.04 


Rivet shear displacement (inch) 


Figure 4.24: Validation example: shear stiffness and strength of rivet spring. 
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0 0.001 0.002 


shear displacement (inch) 

Figure 4.25: Validation example: shear stiffness and strength of adhesive spring. 


global and local modeling levels. The global shell model captures the overall non- 
linear response of the stiffened, curved, pressurized structure. The local shell model 
provides the detailed deformation and stress field near the crack tips to compute 
the fracture parameters (e.g., CTO A) that control stable crack growth. 

4.2.3 Determination of CTOA c 

Flat panel tests were conducted by the Boeing Commercial Airplane Group to 
obtain material properties for fatigue and fracture analysis of the curved fuselage 
panels. Four, 48 inch wide, 80 inch long, 0.063 inch thick middle crack tension 
(MT) specimens were tested. The flat panel specimens were made from the same 
aluminum sheet used for the skin of the curved fuselage panels. A constant ampli- 
tude cyclic loading was applied to propagate an initial sawcut. After the fatigue 
crack growth, a residual strength test was conducted under a monotonically in- 
creasing load. The test matrix prior to the residual strength test is summarized 
in Table 4.1. Visual crack extension measurements were taken. Surface CTOA c 
was measured for Specimen 2024_FAA_TL3 during the residual strength test. Nine 
values were obtained and the mean of the measured critical angles was about 5.5 
degrees with a scatter band about ±1.0°. 

The value of CTOA c used in the residual strength analysis of the fuselage panels 
was determined by finding an angle within the scatter band of the CTOA c mea- 
surements that best correlates with the observed stable crack growth and residual 
strength of the coupon tests. The FRANC3D/STAGS program was used to sim- 
ulate fracture behavior of the MT specimens. A finite element mesh modeling a 
quarter of the specimen with a crack tip element size of 0.04 inch and a half plane 
strain core height equal to 0.08 inch is shown in Figure 4.26. The plane strain core 
was used to capture the three-dimensional (3D) constraint effects developed at the 
local crack tip [94, 31, 56]. The half core height was about the thickness of the 
specimen. 




58 


Table 4.1: Test Matrix for MT Specimens (after [49]) 


Specimen ID 

half initial crack 

half final fatigue crack 

® fatigue 

R 

(inch) 

(inch) 

(ksi) 


2024_FAA_TL3 

2.0 

8.0 

8.0 

0.1 

2024_FAA_TL4 

2.0 

5.5 

16.0 

0.1 


5.5 

8.0 

8.0 


2024_FAA_TL5 

5.0 

12.0 

12.0 

0.1 

2024_FAA_TL6 

2.0 

8.0 

7.0 

0.5 




Figure 4.26: Finite element mesh for a quarter of 48 inch wide MT specimen. 
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o 2024 FAA TL3 

a 2024 I -A A TL4 

predicted results 

.5° 


0 0.5 1 1.5 2 2.5 3 

Half Crack Extension (inch) 

Figure 4.27: Predicted applied stress versus stable crack growth for 48 inch wide 
MT specimen (half plane strain core height = 0.08 inch). 



Figure 4.27 compares the predicted stable crack growth results to the exper- 
imental measurements. The CTOA c of 4.5 degrees best correlates the predicted 
and measured residual strengths. However, it under-estimates the applied stress 
at the early stage of stable crack growth. The 5 and 5.5 degree critical angles give 
a better correlation for the early stable crack growth, but over-predict the residual 
strength by 8.5% and 14.3%, respectively. 

The effect of the 3D constraint zone, i.e., the height of the plane strain core, on 
stable crack growth and residual strength prediction was further investigated. Note 
that in general, a change of the plane strain core height requires a different value 
of CTOA c to correlate the predicted and measured residual strengths. Figure 4.28 
shows the core height effects on the stable crack growth prediction. A slightly 
better correlation for the early growth is observed by increasing the plane strain 
core height. 

The discrepancy between predicted and measured stable crack growth at the 
early stage of tearing might relate to the residual plastic deformation left by the 
fatigue crack growth. This effectively increases the crack opening resistance during 
early stable crack growth [30]. The plastic wake effect on stable crack growth and 
residual strength analysis is further discussed in Section 4.2.5. 

4.2.4 Numerical Results: Comparison with Strain Gage 

Strain gage comparisons were made to verify predicted stress distributions. The 
strain gage readings were recorded during fatigue and residual strength tests. The 
records as the panels reach their residual strengths are of primary interest in this 
study. The corresponding damage configurations are: 
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o 2024 FAA TL3 

a 2024 FAA II A 

— - 5.0° half core height = 0.16 in. 
4.5° half core height = 0.08 in. 

— 4.0° half core height = 0.04 in. 


0 0.5 1 1.5 2 2.5 3 

Half Crack Extension (inch) 

Figure 4.28: Effect of plane strain core height on predicted stable crack growth for 
48 inch wide MT specimen. 

1. a 38.2 inch long crack with a severed central frame for the panel without 
MSD cracks, and 

2. a 41.7 inch long crack with a severed central frame for the panel with MSD 
cracks. 

The ultimate pressure loadings are 9.4 psi and 7.5 psi, respectively. The loca- 
tions of strain gages and the damage configurations are illustrated in Figure 4.29. 
Five back-to-back strain gage rosettes are numbered for the purpose of compari- 
son. Because similar trends for stress distributions were observed for both damage 
configurations, only detailed strain gage comparisons for the panel without MSD 
cracks are described below. 

Nonlinear stress analyses at the global and local modeling levels were per- 
formed. Figure 4.30 shows the overall deformed structures at both levels. Con- 
vergence studies were conducted to ensure accuracy of deformations and stress 
distributions. Figure 4.31 shows three finite element discretizations, Gl, G2, and 
G3, used at the global modeling level. The mesh density around the gage loca- 
tions was progressively refined from global model Gl to G3. The predicted hoop 
stress distributions compared to strain gage readings are shown in Figures 4.32 
and 4.33; the predicted results converge quickly. The predicted membrane hoop 
stresses agreed well with experimental measurements. The predicted bending hoop 
stresses were comparable to experimental measurements as one refined the finite 
element meshes. 

Two discretizations, LI and L2, were performed at the local modeling level 
(Figure 4.34). The mesh density in the local model LI was about the same as the 
corresponding region in the global model G3. The purpose is to ensure transition 






SECTION A-A 


Figure 4.29: Strain gage locations on the skin and tear strap. 
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Figure 4.30: Deformed structures of the validation example at global and local 
modeling levels (pressure = 9.4 psi, crack length = 38.2 inch, magni- 
fication factor = 5.0). 
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inner skin 


outer skin 



Model G1 
Model G2 
Model G3 



Stringer Location 



Stringer Location 



Stringer Location Stringer Location 


Figure 4.32: Global convergence study: comparison between computed and mea- 
sured hoop stresses for strain gage 1-4 (pressure = 9.4 psi; crack length 
= 38.2 in.; frame cut; No MSD). 
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O strain gage measurement 

Model G1 

Model G2 
Model G3 




Frame Location Frame Location 





Frame Location 


Frame Location 


Figure 4.33: Global convergence study: comparison between computed and mea- 
sured hoop stresses for strain gage 5 (pressure = 9.4 psi; crack length 
= 38.2 in.; frame cut; No MSD). 
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Figure 4.34: Two mesh discretizations at the local modeling level. 


accuracy of the hierarchical modeling. Predicted hoop stress distributions from 
local model LI agreed well with global model predictions and experimental mea- 
surements (Figures 4.35 and 4.36). The other discretization, local model L2, had a 
much higher mesh density, which is suitable for stable crack growth analyses. The 
results from local model L2 disagreed with the rest of numerical predictions and 
experimental measurements (Figures 4.35 and 4.36). In particular, the predicted 
membrane hoop stresses were much lower than those observed. 

The discrepancy is related to the idealized representation of the two-noded 
spring element for the rivet connection in the finite element model [150, 141], The 
single point connection results in unrealistic distortion of the surrounding shell 
elements. The local distortion causes premature yielding of the shell elements 
and reduces the load transfer from sheet to rivet. This artificial distortion of the 
shell elements is discretization-dependent [141, pp. 318-327]. Refining the mesh 
captures the local artificial distortion better, but makes the comparison to strain 
gage readings worse [150]. 

Two modeling idealizations are proposed to avoid this artificial effect. One is 
to faithfully represent the geometry of the rivets and their interference with the 
sheets. This will considerably increase the required computational resources and 
may not be simple to implement in thin-shell elastic-plastic crack growth analyses. 
The other approach is to generate distributed connections between the two-noded 
spring element and the surrounding shell elements [150]. The load distribution 
can be accomplished by defining rigid links, stiff spring elements, or least-squares 
loading conditions that connect the rivet-spring node to the surrounding shell- 
element nodes. Care must be taken while defining the area in the shell elements 
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Figure 4.35: Global and local model study: comparison between computed and 
measured hoop stresses for strain gage 1-4 (pressure = 9.4 psi; crack 
length = 38.2 in.; frame cut; No MSD). 
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Figure 4.36: Global and local model study: comparison between computed and 
measured hoop stresses for strain gage 5 (pressure = 9.4 psi; crack 
length = 38.2 in.; frame cut; No MSD). 
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Figure 4.37: Illustration of distributed connection that connects fastener node to 
surrounding shell nodes 


over which the rivet load is distributed. The area should be of the order of the rivet 
cross-sectional area, since distributing the load over a larger area may inadvertently 
stiffen the shell elements. 

Figure 4.37 illustrates simulation of the distributed connection using stiff spring 
elements. The stiff spring elements with an order of magnitude stiffer than the 
rivet spring element were used to distribute the rivet load. For a rivet located 
on a prescribed tearing path, it is expected that the rivet stays intact on only 
one side of the crack as the crack propagates through the rivet. Thus, only the 
shell elements on this side of the crack were used to model the distributed rivet 
connection. Figures 4.38 and 4.39 show the predicted hoop stress distributions with 
distributed connection simulations; a much better prediction is observed. The local 
mesh model, taking into consideration the distributed rivet connection, was used 
for stable crack growth and residual strength analyses. 

4.2.5 Numerical Results: Stable Crack Growth and Resid- 
ual Strength Analyses 

Elastic-plastic crack growth and residual strength analyses were conducted using 
the local model. Both 4.5 and 5.5 degree critical angles computed at 0.04 inch 
behind the growing crack tip were used to investigate the sensitivity of CTOA c on 
stable crack growth and residual strength prediction. The 4.5° CTOA c was the 
angle that best correlates the predicted and observed residual strengths of the MT 
tests. The 5.5° angle was the mean from the surface CTOA c measurements in the 
MT tests. The plane strain core height was 0.16 inch along the prescribed tearing 
path. 

Figure 4.40 shows predicted results from the first attempt for stable crack 
growth analyses. The change of the CTOA c from 4.5° to 5.5° increases predicted 
residual strength by about 33% and 22% for the cases without and with MSD 
cracks, respectively. 



Membrane Hoop Stress (ksi) , Hoop Stress (ksi) 
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Figure 4.38: Effects of distributed rivet connection: comparison between computed 
and measured hoop stresses for strain gage 1-4 (pressure = 9.4 psi; 
crack length = 38.2 in.; frame cut; No MSD). 
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(a) (b) 

Figure 4.40: Comparison between the predicted stable crack growth and experi- 
mental measurements: (a) case without MSD, and (b) case with MSD. 
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Although analysis results in Figure 4.40 clearly demonstrate the loss of residual 
strength due to the presence of MSD, all the predicted results (i) under-estimate 
the pressure loading to initiate the stable crack growth, and (ii) over-estimate the 
residual strength. 

The much lower predicted pressure for tearing initiation is mainly caused by 
residual plastic deformation left by the fatigue crack growth. A possible cause for 
the lower residual strengths observed in the test may be related to the occurrence 
of tear strap failure. Both effects are discussed below. 

4. 2. 5.1 Residual Plastic Deformation Effects 

The test panels were subjected to pressure cycling prior to the residual strength 
test. To incorporate the residual plastic deformations due to the cyclic loading, the 
residual strength analyses were re-performed using an elastic-plastic cyclic loading 
simulation suggested by Newman [89]. The procedure consists of the following 
steps: 

step 1 Close an appropriate length of fatigue crack. 

step 2 Load the fuselage model up to the maximum pressure loading conducted 
in fatigue tests. 

step 3 Release the crack tip node and unload the model. 

step 4 Repeat steps 2 and 3 until the crack tip reaches the initial position for 
stable tearing. 

This procedure implies that the fatigue crack only propagates at the maximum 
pressure during the cyclic loading simulation. For Mode-I only deformations under 
constant-amplitude load cycling, crack surfaces close at a positive applied load (i.e., 
step 3). The contact stresses cause the material to yield in compression. Crack 
face contact and compressive yielding were not modeled in the current simulations. 

In subsequent analyses, the fuselage model is brought to the operating pressure 
level during fatigue tests without allowing the crack to advance. The crack is then 
allowed to advance one element, and the load is returned to zero. Figure 4.41 
illustrates results for a 0.32 inch length of fatigue crack closure used in the analysis 
for the case without MSD cracks. The crack-opening and crack-closure pressures in 
the fuselage panel simulations follow similar trends observed in the MT flat panel 
simulations [89]. After two cycles of simulation, the crack-opening and crack- 
closure pressures quickly stabilize to 7.2 psi and 5.3 psi, respectively. 

Figure 4.42 shows two predicted crack opening profiles with and without fatigue 
crack closure effect when the pressure loading reaches 8.6 psi (no growth). The ef- 
fects of residual plastic deformations on the crack opening profile and consequently, 
the CTOA prediction, are clearly observed. 

The 7.2 psi crack-opening pressure shown in Figure 4.41 seems to be too high in 
comparison with 2D plane stress results [89] and laboratory observations [39, 40]. 
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Figure 4.41: Predicted crack-opening and crack-closure pressure under cyclic load- 
ing (cyclic pressure = 8.6 psi, No MSD). 
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(a) (b) 

Figure 4.42: Predicted crack opening profiles of outer skin at initial tearing crack 
tip: (a) case without fatigue crack closure, and (b) case with 0.32 inch 
fatigue crack closure (no stable crack growth, magnification factor = 
2.0). 

This may be due to the lack of modeling of contact conditions when the crack closes. 
That is, the crack faces pass each other so no compressive yielding is developed 
in the unloaded state. The compressive yielding stress will reduce residual tensile 
plastic deformation thus leading to a lower crack-opening pressure [89]. 

Figure 4.43 shows the predicted results for a 0.08 inch length of fatigue crack 
closure used for the case with MSD cracks. During cyclic loading simulation, the 
lead and MSD crack tips are released simultaneously. The crack-opening and crack- 
closure pressures at the second loading cycle for the lead crack are about 4.7 psi and 
3.3 psi, respectively. We note that the length of fatigue crack closure is restrained 
by the le7igth of MSD cracks. Further amount of fatigue crack closure simulation is 
possible, but leads to somewhat ambiguous MSD fatigue crack propagation. The 
results after two cycles of simulation, however, are believed to essentially capture 
the residual plastic deformation effects. This assertion is based on observations 
from the case without MSD cracks (Figure 4.41). 

Figure 4.44 shows predicted stable crack growth incorporating the closure ef- 
fects. Table 4.2 summarizes the predicted and observed starting pressure to initiate 
stable crack growth. The plasticity-induced closure increases the initiation pres- 
sure by about 150% to 210%. The predicted crack initiation loads are within 6% of 
experimental measurements for the cases that incorporate prior plastic residual de- 
formations due to fatigue crack growth. However, the predicted residual strengths 
are still higher than those observed. 

4. 2. 5. 2 Effects of Tear Strap Failure 

A possible cause for the lower residual strengths observed in the test is the oc- 
currence of failure of other structural elements. Figure 4.45 shows the predicted 
effective stress distribution in the outer skin, inner skin, outer tear strap, and inner 
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Figure 4.43: Predicted crack-opening and crack-closure pressure under cyclic load- 
ing (cyclic pressure = 7.0 psi, MSD). 


Table 4.2: Predicted and Observed Pressure Loading for Stable Tearing Initiation 



predicted (psi) 

observed (psi) 

CTOA c = 4.5° 

CTOA c = 5.5° 

No MSD 

2.3 

2.7 

8.3 

No MSD (0.32 inch closure) 

8.3 

8.4 

MSD 

2.5 

2.8 

6.7 

MSD (0.08 inch closure) 

6.3 

6.5 






operating pressure (psi) 
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Figure 4.44: Comparison between predicted stable crack growth with fatigue clo- 
sure effects and experimental measurements: (a) case without MSD, 
and (b) case with MSD. 
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Figure 4.45: Predicted effective stress distribution (pressure = 9.86 psi, da = 0.5 
inch, CTOA c = 5.5°). 

tear strap as the stable crack growth analysis reaches 9.86 psi pressure loading for 
the case without MSD cracks. Net section yielding is clearly shown in the inner 
tear strap. 

The possible breakage of the inner tear strap during the residual strength test 
was also reported in [49]. To further investigate this damage scenario, a tear 
strap with rivet holes was modeled. By taking the kinematic boundary conditions 
from the local fuselage model, a stress concentration around the holes is observed 
(Figure 4.46). It is then postulated that the high stress concentration is likely to 
initiate new cracks from the rivet holes thus leading to breakage of the inner tear 
strap. 

To incorporate the tear strap damage scenario into the crack growth analysis, 
the inner tear strap is cut prior to fatigue crack closure and stable crack growth 
analyses as illustrated in Figure 4.47. The predicted crack-opening pressures of 
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Figure 4.46: Predicted effective stress distribution of inner tear strap with rivet 
holes (pressure = 9.86 psi, da = 0.5 inch, CTOA c = 5.5°). 

the broken tear strap models with 0.32 and 0.08 inch fatigue crack closure are 7.0 
psi and 3.1 psi for the cases without and with MSD cracks, respectively ( cf . 7.2 psi 
and 4.7 psi for the models with the intact tear strap). 

Figure 4.48 shows the predicted stable crack growth and residual strength for 
the fuselage models with a broken inner tear strap. The predicted residual strength 
using 4.5° CTOA c is within 13% of the experimental observation for the case 
without MSD cracks and within 1% of the experimental observation for the case 
with MSD cracks. 

We further examine several damage scenarios with the possible occurrence of 
the tear strap failure at various stages of stable crack growth. In subsequent 
analyses, the inner tear strap stays intact until it reaches a certain amount of stable 
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Figure 4.47: Illustration of broken inner tear strap. 
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Figure 4.48: Comparison between predicted stable crack growth with broken tear 
strap and experimental measurements: (a) case without MSD, and 
(b) case with MSD. 
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Figure 4.49: Comparison between predicted results with tear strap failure at vari- 
ous stages during stable crack growth and experimental measurements 
(No MSD, CTOA c = 4.5°). 

crack growth. The analysis is continued until the residual strength is reached. 
Figure 4.49 shows the predicted stable crack growth and residual strength using 
4.5° CTOA c for the fuselage models without MSD cracks. For comparison, the 
predicted results with an intact tear strap shown in Figure 4.44 are also plotted. 
The influence of the tear strap failure on residual strength prediction is again 
observed. The occurrence of the tear strap failure at various stages of stable crack 
growth affects the predicted crack growth resistance. But this scenario has a very 
mild influence on residual strength prediction, as long as there is a sufficient amount 
of tearing before the structure reaches its residual strength. 

4. 2. 5. 3 Discussion 

The difference between predicted and observed residual strengths for the case with- 
out MSD cracks may be due to the simulated excess residual plastic deformation 
prior to tearing. One way to reduce the plastic wake is to grow the crack at one half 
the actual fatigue load. The corresponding crack-opening pressure with 0.32 inch 
of fatigue crack closure for the case without MSD is 3.2 psi. This, in conjunction 
with the tear strap damage scenario and 4.5° CTOA c , predicts 9.34 psi residual 
strength for the case without MSD (Figure 4.50). The result is within 1% of the 
experimental observation. However, the crack tearing now initiates at loads much 
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Figure 4.50: Comparison between predicted stable crack growth and experimental 
measurements (broken tear strap, reduced plastic wake, No MSD). 

lower than those seen in the experiment, indicating that this correlation may only 
be coincidental. 

Another possibility of higher predicted residual strength for the case without 
MSD may be related to the fact that the current model does not faithfully model 
fracturing processes in the vicinity of rivets. In the panel test, the lead crack 
propagated into and re-initiated from a rivet hole as illustrated in Figure 4.51. 
Apparently, neither the CTOA fracture criterion for the lead crack propagation nor 
the idealized distributed rivet representation have sufficient accuracy in capturing 
this phenomenon. Further investigation is needed to quantify its effect on residual 
strength prediction. 

4. 2. 5. 4 Major Observations 

Several observations are made from stable crack growth and residual strength anal- 
yses conducted in this section: 

• For all the scenarios simulated, the loss of residual strength due to the pres- 
ence of small MSD cracks is consistently observed. The reduction in residual 
strength caused by MSD varies from 28% to 47%. 

• The residual strength prediction is sensitive to changes in CTOA c . Altering 
the CTOA c from 4.5° to 5.5° changes the predicted residual strength by 17% 
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Figure 4.51: Illustration of crack propagation near rivet: (a) lead crack approach- 
ing rivet, (b) lead crack growing into rivet hole, and (c) new crack 
initiating out of rivet hole. 


to 33% for all the analyses conducted in the damage configuration without 
MSD cracks. It changes the predicted residual strength by 12% to 22% for 
the case with MSD. 

• The residual plastic deformation or the plastic wake from fatigue crack 
growth has a strong effect on stable crack initiation and a mild effect on 
residual strength prediction. For stable crack growth initiation, it is essen- 
tial to incorporate the plastic wake to accurately predict the starting pressure 
loading. Neglecting plastic wake effect leads to a totally erroneous prediction 
of early stable crack growth. For all the residual strength analyses conducted, 
the plastic wake increases the predicted residual strength by 3% to 9%. 

• The breakage of the inner tear strap, categorized as possible failure of other 
structural elements during stable crack growth, is crucial to residual strength 
prediction. For all the analyses conducted, the occurrence of the broken tear 
strap reduces the predicted residual strength by 24% to 30%. Cutting the 
tear strap prior to or during stable crack growth analysis is a preliminary 
approach to model this damage scenario. A better approach would be to 
incorporate proper mechanics to initiate and propagate the damage directly 
in the crack growth analysis. Also, dynamic effects resulting from the failure 
of the tear strap could be simulated, and may not be negligible [130]. 


4.3 Summary 

The crack tip opening angle (CTOA) fracture criterion obtained from coupon tests 
is used to predict fracture behavior and residual strength of built-up aircraft fuse- 
lages that are subjected to widespread fatigue damage (WFD). In the process, 
the feasibility and validity of the analysis methodology are assessed. The major 
findings of the fuselage panel study are: 

1. The occurrence of small MSD cracks substantially reduces the residual strength 
of pressurized fuselages. 
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2. Modeling fatigue crack closure is essential to capture the fracture behavior 
during early stable crack growth. 

3. Possible damage of other structural elements during stable crack growth, e.g ., 
tear strap failure, substantially reduces the residual strength of pressurized 
fuselages. 

4. The distributed rivet load treatment of fasteners is crucial for the local crack 
growth model to accurately predict the stress distribution. 

5. It is apparent that more full-scale fuselage tests need to be conducted to 
fully verify the analysis methodology. However, the CTOA fracture criterion 
together with the FRANC3D/STAGS program proves to be an effective tool 
to simulate: (1) lead crack growth, (2) MSD crack growth, (3) multiple crack 
interaction, (4) plastic wake from fatigue crack growth, (5) tear strap failure, 
and (6) corrosion damage in pressurized fuselages. 



Chapter 5 

Theory for Curvilinear Crack 
Growth in Planar and Thin Shell 
Structures 


Theories and simulations presented in the previous chapters mainly deal with self- 
similar elastic-plastic crack growth where crack trajectories are known a priori. 
However, a crack in a shell-like structural component under combined loading will 
likely propagate in a non-self-similar fashion. Curvilinear crack growth can lead 
to the so-called flapping phenomenon observed in pressurized fuselages [137, 82] 
(Figure 5.1). Flapping can produce a controlled opening in the fuselage that causes 
a “safe” decompression and can prevent catastrophic failure of the structure. 

As discussed in the previous chapters, stress intensity factors (SIFs) and the 
crack tip opening angle (CTOA) serve well to explain fatigue and elastic-plastic 
crack advancement in thin-sheet metallic structures. In addition to these crack 
growth criteria, a criterion for predicting the direction of propagation is needed to 
simulate curvilinear crack growth. 

This chapter together with next two chapters examines some relevant issues 
about curvilinear crack growth simulations. Theories for curvilinear crack growth 
in planar and thin shell structures are discussed in this chapter. In Chapter 6, 
computational methods used to evaluate the T-stress term that is known to have 
a significant effect on crack trajectory prediction are discussed. Curvilinear crack 
growth simulations of coupon tests and full-scale fuselage panel tests are presented 
in Chapter 7. 


5.1 Introduction 

In general, a crack in planar and thin shell structures under mixed-mode loading 
will propagate in a curved fashion. This so-called non-self-similar crack growth 
where crack trajectories are not known a priori requires a direction criterion to 
predict the impending angle of crack propagation. For crack growth in ideally 
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Figure 5.1: Flapping phenomenon in pressurized fuselage due to curvilinear crack 
growth (after [82]). 

brittle isotropic material, the three most common theories for predicting the crack 
propagation angle are: 

1. the maximum tangential stress theory {a eernax ) [41], 

2. the maximum energy release rate theory ( G{6 ) max ) [58], and 

3. the minimum strain energy density theory ( S(9) m i n ) [126], 

There is no consensus on the most satisfactory theory to predict crack growth 
direction. A convenient way to compare the predicted crack growth angles from 
various mixed-mode theories is through introducing a mixed-mode parameter, T: 

T = Ttan -1 | | (5.1) 

7T K n 

that characterizes the elastic loading mixity [124], Comparisons of the elastic 
mixity parameter, T, versus the predicted crack growth angle, 9 C , from the three 
mixed-mode theories are plotted in Figure 5.2. For small values of T where K n is 
dominant, the three theories predict different crack propagation angles. For large 
values of T where Kj is dominant, the three theories predict similar results. In the 
present study, the maximum tangential stress theory is used as a starting point to 
evaluate the direction of crack growth. 

The theory in its original form used only the singular stress fields near the crack 
tip to evaluate the maximum tangential stress [41], Subsequent studies [145, 46] 
suggested that the non-singular stress fields can have a significant effect on crack 
growth direction and crack path stability. Recent studies [152, 72, 103] further 
indicated that the non-singular stress fields play an important role in predicting 
crack turning and flapping in fuselage structures. 
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(Mode II) (Mode I) 

Figure 5.2: Comparison of elastic mixity parameters versus predicted crack growth 
angles. 

The maximum tangential stress theory was originally proposed to predict the 
direction of crack propagation in ideally isotropic brittle material under plane 
stress or plane strain conditions [41]. It has been extended to include the effects 
of elastic and fracture anisotropy [14, 12, 21, 104], Boone et al. [12], using the 
theory proposed by Buczek and Herakovich [14], showed that both elastic and 
fracture orthotropy can affect the direction of crack propagation, but the fracture 
orthotropy was found to be a much more dominant factor. Chen et al. [21] found 
that the fracture orthotropy has a strong influence on predicted crack trajectories 
in narrow body fuselages. The theory has recently been extended to include the 
non-singular stress contributions [104], 

All the above theories were developed for two-dimensional, linear elastic frac- 
ture mechanics (LEFM) problems. As discussed in Section 1.2.2, for pressurized 
thin shell structures a geometrically nonlinear analysis is required to capture the 
crack tip deformations. The fracture parameters developed under the linear elastic 
framework have been extended to handle geometrically nonlinear problems with 
finite elastic deformations. Eshelby [43], using the energy-momentum tensor, de- 
veloped a Lagrangian framework for LEFM problems. The counterpart fracture 
parameters in the Lagrangian formulation are able to characterize the crack tip 
fields for deformations of arbitrary magnitude [43, 119, 125]. 

The maximum tangential stress theory originally developed under small-scale 
yielding conditions has been extended to the elastic-plastic range. Shih [124] stud- 
ied mixed-mode, plane strain problems using the deformation plasticity and near- 
field singularity dominated by Hutchinson-Rice-Rosengren (HRR) fields [60, 120, 
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Figure 5.3: Local crack tip coordinates. 


59]. He concluded that the crack propagation angle depends not only on the elas- 
tic mixity parameter, T, but also on the strain hardening exponent of materials. 
Maccagno and Knott [80] as well as Pawliska et al. [100] using Shih’s approach, 
further showed the applicability of the maximum tangential stress criterion for 
elastic-plastic materials. Recent studies [134, 51, 133] indicated that an additional 
directional criterion related to shear type fracture is needed for Mode-II dominated 
crack propagation in metals. These authors, however, agreed that the maximum 
tangential stress criterion suffices for Mode-I dominated crack propagation. 

5.2 Crack Tip Fields in Two Dimensions and Thin 
Plates 

The crack tip stress and displacement fields in two dimensions [146] as well as in 
thin plates subjected to bending and twisting [147, 57] are outlined below. 

Let ( x , y ) be the local Cartesian coordinates and (r, 9) be the local polar co- 
ordinates centered at the crack tip (Figure 5.3). For two-dimensional elastic crack 
problems, Williams [146] derived a set of solutions for stresses and displacements 
that would satisfy equilibrium and compatibility equations in the neighborhood of 
a crack tip: 


+oo 


°ij = Axr * fijW 

(5.2) 

A=— oo 


+oo 


u i = Bx r * +1 9i(°) 

(5.3) 


A=— oo 


where | is the eigenvalue of the problem and .4 A and B\ are coefficients of expan- 
sions. 

With the physical argument that the total strain energy should be bounded at 
the crack tip, stress expansions from Equation (5.2) in terms of the local Cartesian 
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where Ki and K n are Mode I and Mode II stress intensity factors, T is the constant 
non-singular stress term that appears only in the a xx component, and A Jn and A IIn 
are Mode I and Mode II coefficients of higher order terms. 

Similarly, displacements from Equation (5.3) can be expressed as: 
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where E is the Young’s modulus, G is the elastic shear modulus, and v is the 
Poisson ratio. For plane stress problems, k = (3 — z/)/(l + i/), E = E, and u = v. 
For plane strain problems, k = (3 — 4i/), E = E/( 1 — i/ 2 ), and v = vj{ 1 — u). 

In addition to the above two-dimensional fields, Williams [147] and Hui and 
Zehnder [57] further derived an asymptotic field for bending in elastic thin plates. 
The bending stress and displacement fields near the crack tip in the context of 
Kirchhoff plate theory are: 
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Figure 5.4: Local crack tip coordinates for a through crack in a plate (after [57]). 
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(5.12) 


where (r, 9, z) are the polar cylindrical coordinates centered at the crack tip 
(Figure 5.4), t is the plate thickness, and ki, k? are the Ivirchhoff theory stress 
intensity factors. 

For elastic curvilinear crack growth in thin shell structures, the local crack tip 
fields are assumed to be sufficiently characterized by the two-dimensional plane 
stress and the Kirchhoff plate fields [57, 108]. For points on a shell midsurface (i.e., 
z = 0), k\ and ko, stress intensity factors make no contributions to the displacement 
and stress fields. We thus will directly extend the two-dimensional crack growth 
direction criterion to handle thin shell LEFM problems. 


5.3 Crack Growth Direction Criterion Based on 
Maximum Tangential Stress Theory 

In their work on predicting the direction of crack growth, Erdogan and Sih stated 
[41]: 

“ . . . (There are) two commonly recognized hypotheses for the exten- 
sion of cracks in a brittle material under slowly applied plane loads: 

(a) The crack extension starts at its tip in radial direction. 

(b) The crack extension starts in the plane perpendicular to the direc- 
tion of greatest tension. 

These hypotheses imply that the crack will start to grow from the tip 
in the direction along which the tangential stress ago, is maximum . . . 


The tangential stress agg near the crack tip can be derived from the local Carte- 
sian stresses (Equations (5.4), (5.5), and (5.6)) with coordinate transformations. 
For two dimensional mixed-mode problems, agg up to the order of the T-stress 
term is: 

ago = . cos - ( K[ cos 2 - — -Ku sin 9 ] + T sin 2 9 (5.13) 

v27rr 2 \ 2 2 J 
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Based on the hypotheses that crack extension would take place in the direction 
along which <j ee possesses a maximum value, we have: 


does n , d 2 aeo 

~af =0 and »r <0 

Taking the derivative of oee with respect to 6 , we have: 
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Taking the second derivative, we have: 
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Thus, gqq is a relative maximum if the following conditions are satisfied: 


{ 



(— | cos (Kj sin0 c + K ir (3 cos 9 C — 1)) + 2 T sin 9 C cos 9 C = 0 
(— | (3 cos 9 C — 1)) (A'/ cos — Kji sin + 2 T cos 2 9 C < 0 


(5.17) 


where r c is a critical distance away from the crack tip and 9 C is the corresponding- 
crack propagation angle. We note that this directional criterion is the same as 
the one proposed by Williams and Ewing [145] and later corrected by Finnie and 
Saith [46]. The above equation reduces to the classical Erdogan and Sill directional 
criterion [41] if T = 0 or r c = 0, i.e., 


Ki sin 9 C + Kn (3 cos 9 C — 1) = 0 


(5.18) 


By comparing Equation (5.17) with Equation (5.18), one finds that a length 
parameter, r c , is needed to incorporate T-stress into the crack growth directional 
criterion. We will further discuss the physical meaning of r c in Section 5.3.3. 


5.3.1 Crack Path Stability under Pure Mode I Conditions 

Using Equation (5.17), crack path instability in pure Mode I conditions is predicted 
for certain circumstances under positive T-stress. For pure Mode I problems where 
K n = 0, we have closed-form solutions for Equation (5.17) [152]: 

e c = 0 (5.19) 


or 


9 C = 2 cos 


-l 


S F ± \l - + Sp\ 


(5.20) 
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Figure 5.5: The tangential stress o ee distributions for different values of critical 
distance r c under pure Mode I conditions (T > 0). 


in which Sf = (3 K r )/ (32T v / 27rr c ). For 9 C = 0, the second derivative of gqq with 
respect to 9 is: 


0 2 o ee 
09 2 


3 K, 

4 y/2l vF c 


+ 2 T 


(5.21) 


Thus, under a negative T-stress environment, 0 2 ( 700 / 09 2 is always smaller than 
zero. This implies that 000 at 9 C = 0 is a relative maximum and the crack will 
grow in a straight (i.e., self-similar) manner. For positive T-stress, the crack will 
propagate in a self-similar fashion only if: 


128tt T J 


(5.22) 


For r c > (9A'|)/(1287tT 2 ), a 00 at 9 C = 0 is a relative minimum and a crack will 
grow in the direction predicted by Equation (5.20). Figure 5.5 illustrates the 000 
distribution for different values of r c and the relative maximum and minimum. 


5.3.2 Determine Crack Propagation Angle under General 
Mixed-Mode Problems 

For general two-dimensional mixed-mode problems, the crack propagation angle 
can be solved from Equation (5.17) for given r c , K I: K n , and T. Figure 5.6 plots 
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T = Hy/2^ 


Figure 5.6: Predicted propagation angle from maximum tangential stress theory 
(Equation (5.17)). 

the predicted propagation angle 0 C versus a dimensionless parameter T [104]: 

T = ^^ry/2nr c (5.23) 

for various ratios of K n /Kj. 

The parameter r c is assumed to be a material constant and will be further 
discussed in next section. A'/ and Kjj for two-dimensional problems can be accu- 
rately computed from the path independent M-integral [149] or from symmetric 
and anti-symmetric parts of the J-integral [15]. For thin shell problems, mem- 
brane and bending stress intensity factors can be obtained from an extension of 
the modified crack closure integral method [106, 142], How to obtain an accurate 
T-stress numerically is not obvious from the literature. We will further discuss this 
issue in Chapter 6. 

It is of historical interest to consider a specific angled crack problem shown in 
Figure 5.7. For this special case, we have [41, 46]: 

Ki = sin 2 (3 K n = sin (3 cos [3 T = cr(cos 2 (3 — sin 2 (3) 

The predicted propagation angle Q c can be solved by applying Equation (5.17). In 
Figure 5.7, the predicted propagation angles 9 C from Equation (5.18) (or from Er- 
dogan and Sill [41]) are compared with those of \Jlr c ja = 0.1 from Equation (5.17) 
(or from Finnie and Saith [46]). Some results are observed: 
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Figure 5.7: Predicted propagation angle for the angled crack problem. 

• Two curves intersect at /3 = 45° in which T = 0. 

• The incorporation of the T-stress predicts a larger propagation angle under 
positive T-stress conditions. It predicts a smaller angle under negative T- 
stress conditions. 


5.3.3 Discussion: The Critical Distance r c 

As shown in Figures 5.5 and 5.6, the critical distance ahead of the crack tip, 
r c , plays an important role in predicting crack path stability and crack growth 
direction. Earlier studies by Williams and Ewing [145] and Finnie and Saith [46] 
assumed that rja is a constant value (cf. 2 r c /a = 0.1 in Figure 5.7). Streit 
and Finnie [132] further postulated that r c is a critical distance in front of the 
crack tip where fracture would occur. A photoelastic and experimental study 
of Mode I crack extension was conducted. The ratio Ki/T was determined by 
analyzing the isochromatic-fringe geometry in the photoelastic experiments. The 
critical distance r c was determined by simply observing the onset of crack turning 
where r c = (9A'|)/(1287tT 2 ) at this instance. They concluded that for 7075-T651 
aluminum plate, r c = 0.005 inch for side-grooved specimens and r c = 0.01 inch for 
the ungrooved, L-T specimens. These values seem to be too small in comparing 
with subsequent experimental studies [110, 103]. 

Ramulu and Kobayashi [110] extended the method of Streit and Finnie [132] 
and measured r c for a dynamically growing crack. They observed that r c is a 
constant value for specimens under various mixed-mode conditions. Based on the 
dynamic photoelastic experiments, they concluded that r c is about 0.05 inch for 
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Homalite-100. Recent studies by Kobayashi and his associates [71, 72, 73] further 
estimated that r c is about 0.04-0.06 inch for 7075-T6 and 2024-T3 thin-sheet 
fracture specimens. 

Pettit et al. [103] determined the magnitude of r c by analyzing crack turning 
in 2024-T3, double cantilever beam (DCB) specimens. They estimated that the 
value of r c is about 0.05 inch for stable tearing under a monotonically increasing 
load. For slow fatigue crack growth under a low stress level of cyclic loading, r c is 
assumed to be negligible. Recent ongoing research [102] further suggested that the 
magnitude of r c for stable tearing in 2024-T3 specimens is about 0.04-0.09 inch. 

In this study, we simply assume that r c is a material constant that can be deter- 
mined from fracture coupon tests. The effect of r c on crack trajectory prediction 
will be extensively examined in Chapter 7. 


5.3.4 Anisotropic Problems 


Material grain structure variations and other micro-structural factors from differ- 
ent forms of material processing can influence fracture toughness variation with 
direction and, therefore, the crack growth trajectory. Taking a rolled sheet made 
of 2024-T3 aluminum alloy for example, the direction perpendicular to the rolling 
direction could have a 5 to 20% higher toughness than that of the rolling direction 
[47, 44, 99]. 

In general, crack propagation in anisotropic media is considerably more com- 
plicated than the isotropic case [12]. In the present work, a simple extrapolation 
of the maximum tangential stress theory to materials with orthotropic toughness 
proposed by Buczek and Herakovich [14] is used. The tangential stress is normal- 
ized with respect to the directional strength of the material. Crack propagation is 
assumed to be in the direction of maximum normalized stress, such that: 


Maximum 


voeiki, k-n, T,r c , 9) 
K c (a) 


ojw 

Kr 


at 9 = 8 r 


critical 


(5.24) 


where a is the angle characterizing the material grain orientation, K c (a ) is the 
strength parameter characterizing the material fracture resistance, and 9 C is the 
angle of impending crack propagation. 

A simple elliptical function is used in this study to characterize the anisotropic 
fracture toughness, K c (a ) [14, 68]. The equation of the ellipse with fracture tough- 
ness K c ( 0°) along the material longitudinal (x) direction and A' c (90°) along the 
transverse ( y ) direction is (Figure 5.8): 


[A- c (0«)2] + [ii- c (90") 2 ] 


(5.25) 


Substituting x 


K c {a) cos cr and y = K c (a ) sin a into Equation (5.25), we have: 


I<c(af 


COS 2 O' 


+ 


sin 2 a 


K c ( 0°) 2 K c (90°)' 


= 1 


(5.26) 
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Figure 5.8: Elliptical representation of anisotropic fracture toughness. 
Thus, the anisotropic fracture toughness K c (a) can be expressed as [104]: 


K c (a) 


cos 2 a T 


1 

d/A',,,)- sin 1 ’ o 


(5.27) 


where K m is the fracture orthotropy ratio dehned as K rn = K c (90°)/ K c (0°). 

Both the fracture orthotropy ratio, K m , and the material orientation angle, a, 
can affect the predicted angle of impending fracture propagation, 9 C , as demon- 
strated in Figure 5.9 and 5.10. 

We note that Equation (5.26) can be generalized to the n-th order: 


and n 


( cos 2 a: i sin 2 a ^ 

\K c {0°) n + K c {90°) n ) 


— 1 was used in [14, 12]. 


1 


(5.28) 


5.4 Discussion: Crack Growth Direction Crite- 
rion for Geometrically and Materially Non- 
linear Problems 

Possible extensions of the above crack growth directional criteria developed un- 
der the LEFM framework to handle geometrically and materially nonlinear shell 
problems are discussed below. 
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5.4.1 Geometrically Nonlinear Problems 


To compute stress intensity factors and T-stress with elastic deformations of ar- 
bitrary magnitude, one can simply evaluate the fracture parameters based on the 
computed quantities in the current deformed equilibrium configuration. This so- 
called Eulerian approach together with the modified crack closure integral method 
has been used successfully to compute the membrane and bending stress intensity 
factors for geometrically nonlinear shell problems [106, 142]. 

It is often desirable to use the path independent integral to evaluate the fracture 
parameters. However, the integral is difficult to be directly applied in the Eulerian 
framework [43]. Taking the deformed structure shown in Figure 4.30 for example, 
it is conceptually difficult to evaluate a rigorous path independent integral at the 
deformed configuration with the occurrence of severe bulging. 

An alternative way is to evaluate the path independent integral in the La- 
grangian framework [43, 70]. The derivations rely on finding the counterparts of 
conservative (i. e., path independent) integrals, well-defined under elastic states 
with infinitesimal deformations, in the context of finite elastic deformations. The 
fracture parameters (for example, stress intensity factors) are then related to these 
conservative integrals. 

The well known conservative J-integral in two dimensions, for example, is given 
by: 

J =1 

r 


' jrr c 9Ui ' 

H o X j on 




n.j d, 


(5.29) 


where , is an arbitrary counter-clockwise contour around the tip of a crack, W 
is the strain energy density, 5 X j is the Kronecker delta, are components of the 
Cauchy stress tensor, tq are components of the displacement vector, and rij are 
components of the normal vector along the contour ,. 

The counterpart of the J-integral under finite elastic deformations can be de- 
rived in a relatively straightforward manner in Lagrangian coordinates, X. With 
re-interpretation of the field quantities with reference to the undeformed config- 
uration, the J-integral for geometrically nonlinear problems can be expressed as 
[43, 70]: 

Jgn = 

c 


wx dUi 

WSxj Pij g 


NjdC 


(5.30) 


where C and Nj are evaluated in the undeformed configuration, W is interpreted 
as the strain energy per unit undeformed volume, and p t j are components of the 
nominal stress tensor (transpose of the first Piola-Kirchhoff stress tensor). 

The J in linearized and Jgn in finite elastic states both characterize the energy 
release per unit crack advance [43]. Under Mode-I, plane stress conditions, we thus 
have: 


(5.31) 


K, = VEJ 
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and equivalently, 

{Ki)gn = \J EJ gn (5.32) 

Similarly, other fracture parameters can be related to their finite deformation 
counterparts. As a result, crack propagation and direction criteria derived under 
LEFM can be extended to handle geometrically nonlinear problems. 

However, we note that the Lagrangian counterparts of conservative integrals 
to characterize the stress intensity factors and T -stress for geometrically nonlinear 
shell problems are yet to be derived. In this study, we will simply use the mod- 
ified crack closure integral method to compute the membrane and bending stress 
intensity factors [106, 142] and the displacement correlation method to evaluate 
the T-stress for geometrically nonlinear shell problems [69, 134, 104], 

5.4.2 Materially Nonlinear Problems 

The crack growth directional criteria developed above are strictly valid only for 
small-scale yielding problems. The criteria may be sufficient for curvilinear crack 
growth under fatigue loading (c/. Section 1.2.2). For stable crack growth, direc- 
tional criteria need to be extended to the elastic-plastic range. 

In an early attempt, Shill [124], using HRR fields, extended the maximum 
tangential stress theory to the elastic-plastic range under plane strain conditions. 
Figure 5.11 shows the predicted crack propagation angle, 9 C , versus the elastic load 
mixity, T (Equation (5.1)). The results indicate that the crack propagation angle, 
0 C , depends not only on the elastic mixity parameter, \k, but also on the strain 
hardening exponent of materials, n. 

For n = 1 (i. e., linear elastic material), the criterion reduces to the maximum 
tangential stress theory of Erdogan and Sih [41], Thus, the theory will not capture 
the crack path instability under the pure Mode-I conditions caused by T-stress. 
And more importantly, the theory is based on the HRR fields. The fields are 
no longer valid to characterize the elastic-plastic crack tip fields after a sufficient 
amount of stable crack growth. These two main disadvantages prohibit direct 
application of the above theory to predict the direction of stable crack growth in 
fuselage structures. 

Recently, a simple crack growth directional criterion based on the crack tip 
opening displacement (CTOD) concept has been proposed by Sutton et al. [133]. 
The criterion is motivated by the laboratory observations of recent Arcan specimen 
tests conducted by Amstutz et al. [1, 32], The test results show that there is a 
sharp transition of crack growth behavior from predominantly Mode I type to 
Mode II type fracture for 2024-T3 thin sheet materials. Since the crack growth 
direction prediction based on the maximum tangential stress theory is mainly for 
Mode I dominated fracture, Sutton et al. [133] proposed a general CTOD-based 
criterion to overcome this disadvantage. 

The CTOD-based crack growth directional criterion considers a kinked crack 
departing from a main crack as shown in Figure 5.12. The criterion postulates 
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Figure 5.11: Comparison of elastic mixity parameters versus predicted crack 
growth angles for different strain hardening exponents (after [124]). 


that crack growth direction of a kinked crack would likely follow a path that gives 
the maximum CTOD of the kinked crack (i.e., 6). The following steps outline 
the procedure to incorporate the directional criterion into the elastic-plastic crack 
advancement controlled by the CTOA criterion [133]: 

step 1 Apply the load monotonically until D (or CTOD) at a specified distance L 
behind the main crack tip, or equivalently CTOA, reaches its critical value. 

step 2 Release the main crack tip node and extend the crack along every possible 
ray from the crack tip. 

step 3 Compute 5 of all possible kinked cracks. 

step 4 Determine the crack propagation angle 9 C by finding the path that gives 
the maximum 5. 

step 5 Continue stable crack growth simulation controlled by the CTOA criterion. 

The CTOD-based directional criterion predicts both Mode I and Mode II type 
crack growth observed in Arcan specimens [133]. The directional criterion, how- 
ever, suffers from its computational inefficiency since 5 must be evaluated for a 
large number of rays to determine the propagation angle. 

To overcome the drawback, Sutton [133] further assumed that there exists a 
unique relationship between S (CTOD of the kinked crack) and D (CTOD of the 
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main crack kinked crack 

Figure 5.12: Illustration of main and kinked crack relationship for CTOD-based 
crack growth directional criterion. 
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Figure 5.13: Comparison of predicted propagation angles for the angled crack prob- 
lem. 


main crack). By defining a local mixity a = D y /D x of the main crack, an empirical 
equation to predict the crack propagation angle, 9 C , was found: 


0 C = 


«i tan 1 (b\Oi) 

a 2 cos(6 2 a)]^| 


if a < O'c 
if a > a c 


(5.33) 


where a c is the critical local mode mixity for the transition between Mode I and 
Mode II type fracture and a i, &i, a 2 , and b 2 are the curve fitting parameters based 
on fracture test data. A set of parameters, a c = 70°, Gq = —36.5, bi = 2.2, 
a 2 = 57.3, and b 2 = 1.0 were obtained from the test data of 0.09 inch thick Arcan 
specimens made of 2024-T3 aluminum. 

Equation (5.33) and the fitting parameters based on Arcan test data were 
argued to be material constants and to be applicable to other geometries with the 
same thickness [133]. We, however, found that this assertion may not hold for all 
cases. Taking the angled crack problems under LEFM for example, a substantial 
difference as shown in Figure 5.13 is observed between the predicted angles by the 
maximum tangential stress theory and those from Equation (5.33). The difference 
observed in Mode I dominated fracture is thought to be related to: (1) built-in 
orthotropy and (2) widespread plasticity in Arcan specimens. Further investigation 
is needed to fully justify the observation and assess the geometry independence of 
Equation (5.33). 

Neither the HRR-type extension of the maximum tangential stress theory nor 
the CTOD-based directional criterion seems to be sufficient to fully characterize 
the direction of elastic-plastic crack growth. In this study, we will simply apply 
the LEFM approach to predict the direction of elastic-plastic crack growth. 
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5.5 Summary 

In this chapter, theories for curvilinear crack growth in planar and thin shell struc- 
tures were discussed. A well known crack growth directional criterion based on the 
maximum tangential stress theory was examined. Singular as well as non-singular 
constant stress (or T-stress) fields were included in evaluating the stress states near 
the crack tip. 

Equation (5.17) described the predicted impending propagation angle based 
on the maximum tangential stress theory with the T-stress effect. Its numerical 
outcome on crack path instability and crack growth direction was discussed in 
Sections 5.3.1 and 5.3.2. This criterion will be used in Chapter 7 to predict the 
propagation angle for isotropic mixed-mode problems. 

In Section 5.3.4, the directional theory was further extended to include the 
effect of material anisotropy. A simple elliptical function was used to characterize 
the anisotropic fracture resistance in different material orientations. The predicted 
propagation angle incorporating the T-stress and fracture toughness orthotropy ef- 
fect was described symbolically in Equation (5.24). The effect of fracture toughness 
orthotropy ratio and the material orientation angle on the predicted propagation 
angle was shown in Figures 5.9 and 5.10. This directional criterion will be used in 
Chapter 7 to predict the propagation angle for orthotropic mixed-mode problems. 

We then discussed possible extensions of the above directional criteria to han- 
dle both geometrically and materially nonlinear problems. For elastic deformations 
of arbitrary magnitude, one can evaluate the stress intensity factors and T-stress 
based on either the Eulerian or Lagrangian formulation. The latter is conceptually 
simple to be used with the powerful path independent integral. However, the La- 
grangian counterparts of conservative integrals to characterize the stress intensity 
factors and T-stress for geometrically nonlinear shells are yet to be derived. In this 
study, the modified crack closure integral method and the displacement correlation 
method will be used to evaluate the stress intensity factors and T-stress in thin 
shells, respectively. 

The possible extension of the maximum tangential stress criterion to predict the 
propagation direction of elastic-plastic crack growth was commented and critiqued. 
A new CTOD-based direction criterion was also examined. We concluded that 
neither of them seems to be sufficient to fully characterize the direction of elastic- 
plastic crack growth. Thus, in this study we will simply apply the LEFM approach 
(Equations (5.17) and (5.24)) to predict the elastic-plastic crack growth direction. 



Chapter 6 

Numerical Evaluation of T-Stress 


In this chapter, numerical methods for T-stress evaluation are discussed. Among 
all the possible methods that can be used to compute T-stress, we focus on the 
path independent integral method because its inherent nature allows us to evaluate 
the desired value in a far-held region away from the crack tip where numerical 
accuracy is greater. We will first put the FRANC3D/STAGS program aside and 
use a powerful two-dimensional hp- version finite element code [74] to fully quantify 
and assess the accuracy of computed T-stress using the path independent integral 
method. We will then discuss applicability of the FRANC3D/STAGS program in 
evaluating T-stress for two-dimensional as well as thin-shell problems. 

6.1 Introduction 

The second term of the elastic asymptotic stress series near a crack tip [146], often 
called T-stress, is known to have significant influence on crack growth direction 
and crack path stability [46, 103, 26]. In addition the T-stress is also known to 
have a strong influence on crack-tip constraint [38, 98]. To obtain an accurate 
T-stress for complex geometries subjected to arbitrary loading thus becomes an 
important task for fracture analysis assessment. 

Several numerical methods have been used to evaluate the T-stress [76, 77, 67, 
123, 129]. An earlier study by Larsson and Carlsson [76] determined the T-term 
from two finite element solutions, one with a A- field and the other with actual load- 
ing and geometry configurations. Leevers and Radon [77] computed the T-stress 
by incorporating the eigenfunctions from Williams [146] in a variational formula- 
tion. Sham [123] used second order weight functions through a work-conjugate 
integral to calculate the T-term. 

To compute T-stress in conjunction with finite element analyses, a path inde- 
pendent integral similar to the J-integral [42, 114] for stress intensity factors is 
highly desirable. Cardew et al. [16] and Kfouri [67] presented a novel J-integral 
type of path independent integral for computing the T-term from finite element 
analyses. Recently, another type of path independent integral for T-stress compu- 
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tation based on the Betti-Rayleigh reciprocal theorem has been proposed [129, 151]. 
Due to its simplicity, the Betti-Ravleigh reciprocal type of conservative integral has 
been widely used to solve crack and notch problems for homogeneous as well as 
bimaterial bodies [131, 127, 8]. As will be shown in this study, these two path 
independent integrals are analytically equivalent. 

To our best knowledge, none of the previous studies has fully addressed the 
accuracy of numerical T-stress computations. Published values vary between three 
and five percent for identical loading and geometry configurations and the error 
for the computed T-stress is generally unknown. R is well known that to increase 
the accuracy of a finite-element computation, either the mesh has to be refined 
(h-version) or the polynomial degrees of the shape functions have to be increased 
(p-version). A combination of both strategies, referred to as the hp- version of 
the finite element method, is known to show exponential rates of convergence in 
energy norm even if the problem has singularities [6]. In our contribution we will 
show that path independent integrals, in conjunction with hierarchical, p- and 
hp-version finite element methods [141], provide a powerful tool to obtain highly 
accurate numerical results for T-stress. Using a novel error estimator for the T- 
stress, the accuracy of the computation is quantified and assessed. 

6.1.1 Outline for Numerical Assessment of T-stress Com- 
putation Using a p-version Finite Element Method 

The derivations of path independent integrals for T-stress are studied in Sec- 
tion 6.2. Finite element implementations of equivalent domain integrals in con- 
junction with the hierarchical p-version finite element method are discussed in 
Section 6.3. To quantify the error in computing T-stress using the path indepen- 
dent integrals, an error estimator is proposed. Using a highly accurate hp-version 
finite element code [74] a benchmark example with various K r and T imposed 
boundary conditions is studied in Section 6.4 with the goal being to assess the 
accuracy of the numerical computation of T-stress. We then compute values of 
T-stress for various well known fracture specimens and compare our results with 
values from the literature. 


6.2 Path Independent Integral For T-Stress Com- 
putation 

Two types of path independent integrals for T-stress evaluation have recently been 
proposed. One is based on the Betti-Rayleigh reciprocal theorem [129, 151] and 
the other is based on Eshelby’s energy momentum tensor [16, 67]. For the Betti- 
Rayleigh reciprocal type of conservative integral, we will detail the derivation be- 
cause of its relative ambiguity in the literature. For Eshelby type integrals, we will 
briefly outline the formulation in Cardew et al. [16] and discuss their analytical 
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equivalence to the Betti-Rayleigh reciprocal type of conservative integral. 

6.2.1 Construction of Path Independent Integral for T- 
Stress Using Betti-Rayleigh Reciprocal Theorem 

The Betti-Rayleigh reciprocal theorem states that “for a hyperelastic body subject 
to two infinitesimal systems of body and surface forces, the work done by the first 
system in the displacement caused by the second equals the work done by the 
second in the displacement caused by the first” [83]. For elastostatic problems, 
the body forces fi and surface tractions U produce displacements tq. The body 
forces f* and surface tractions t* produce displacements u* and are called auxiliary 
fields. From the results of divergence theorem, assuming sufficient smoothness of 
the functions and the boundaries of the body, we can prove the Betti-Rayleigh 
reciprocal theorem: 

t* Ui dA + HI ^ f* Ui dV =11 t iU * dA + HI fin* dV (6.1) 

where A is the bounding surface of the body, dA is an infinitesimal element of the 
surface, V is the volume of the body, and dV is an infinitesimal element of the 
volume. For a two-dimensional case without body forces, Equation (6.1) reduces 
to: 

- Uu*) dS = 0 (6.2) 

s 

where S is the bounding curve of the body and dS is an infinitesimal segment of 
the curve. 

To construct a path independent integral for a two-dimensional elastic body 
with a crack using the Betti-Rayleigh reciprocal theorem, the procedure outlined 
by Stern et al. [131] is followed. First consider a contour integral along a closed 
path (C, C+, C e , and C_) as shown in Figure 6.1. From Equation (6.2), we have: 



tjUj )d,C + + 



tiU*)dC e 


+ 



tiU*)dC _ 


0 


(6.3) 


Since C + and C- are traction free, we have 



tiU*)dC e 


(6.4) 


Equation (6.4) proves the path independence of the contour integral. With l t = 
where o rj are components of the stress tensor and n :j are components of the 
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Figure 6.1: A closed contour in the neighborhood of a crack tip. 


outward pointing normal vector of the contour, Equation (6.4) can be expressed 
as: 



(a-jUi - (Ji j u*)n j dC 



(a-jiii - a ij u*)n j dC e 


(6.5) 


Due to its path independent nature, the integral on the left hand side of Equa- 
tion (6.5) can be evaluated on a contour away from the crack tip where numerical 
solutions can be used for and u t without losing too much accuracy. The integral 
on the right hand side of Equation (6.5) is evaluated analytically as e— >0 [8]. 

Let ( x , y) be the local Cartesian coordinates and (r, 9) be the local polar 
coordinates centered at the crack tip. For two-dimensional elastic crack problems, 
Williams [146] derived a set of solutions for stresses and displacements that satisfy 
the equilibrium and compatibility equations near a crack tip ( cf . Equations (5.2)- 
(5.8)): 


+oo 

a ij = fiji 9 ) ( 6 ‘ 6 ) 

A=— oo 
+oo 

Ui = Bx r * +1 9*(°) ( 6 - 7 ) 

A=— oo 


where | is the eigenvalue of the problem and A\ and B\ are coefficients of the 
asymptotic expansions. In order to obtain coefficients of a particular order of ~ 
alone, the auxiliary fields required in the reciprocal work relation in Equation (6.5) 
are: 

a*. ~ r~ 2 ~ 2 v* ~ r - ^- 1 (6.8) 

The question arises how to extract the contributions of T-stress from above se- 
ries expansions, without contributions from singular aud other higher order terms. 
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Considering Equation (6.8), the idea is to choose auxiliary fields so that cq* ~ r~ 2 
and u* ~ r _1 to cancel all contributions from the first, singular terms of the 
expansions in Equation (6.5) to the T-stress term. The auxiliary stresses and 
displacements in local Cartesian coordinates are: 


* 

Cx.T 

cos 29 + cos 4 9 
2 ti r 2 

(6.9) 

* 

°yy 

cos 29 — cos 4 9 
2 ti r 2 

(6.10) 

* 

Try 

sin 4 9 

2 ti r 2 

(6.11) 

v* — — 

1 K, cos 9 + cos 3 9 

(6.12) 

a x ~ 

4 ti r 2 G 

v* — - 

1 — « sin# + sin 30 

\/ 

(6.13) 

u y ~ 

4 ti r 2 G 


where G is the shear modulus, and k = ( Z — v)/(l + v) for plane stress problems 
and k = (3 — Av) for plane strain problems. 

It is clear that with such auxiliary fields, other higher order terms in stress and 
displacement expansions result in no contribution to the contour integral as r— >• 0. 
With some algebraic manipulation, it is possible to show that no contribution 
occurs from the singular terms, and 

T = E (cr^Ui — <JijU*)rijdC e as e— >-0 (6.14) 

JCe 

where E = E for plane stress problems and E = E/( 1 — v 2 ) for plane strain 
problems in which E is Young’s modulus. 

As a result, T-stress is readily computable by combining Equations (6.5) and 
(6.14) with finite element analyses: 

T = E {<J*ju[ E - of E u*)rijdC (6.15) 

where uf E and af E are stresses and displacements of a finite element solution. 

6.2.2 Construction of Path Independent Integral for T- 
Stress Using Eshelby’s Energy Momentum Tensor 

Another type of path independent integral, following Eshelby, has been proposed 
for T-stress computations [16, 67]. The formulation uses auxiliary fields from point 
force loading in conjunction with finite element results. The T-stress is obtained 
by combining a common ./-integral with a /-integral of superimposed auxiliary 
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Figure 6.2: Point force applied at crack tip. 


fields from the point force solutions. A brief outline of the formulation is given 
below. 

Considering a point force at the crack tip in a infinite body (Figure 6.2), the 
stress fields in local Cartesian coordinates are: 


O xx 

a yy 

®xy 


cos 3 9 

it r 

f 

— — cos 9 sin 2 9 
nr 


f 

cos 2 9 sin 9 

7 r r 


(6.16) 

(6.17) 

(6.18) 


Let F denote the two-dimensional elastic solution near the crack tip and F the 
solution from the point force. The J-integral of the superimposed state of F and 
F can be expressed as: 


J(F, F) 



ip%k T Cbfc)(f?ifc + €ik)3-xj 


{< 7 jj + + U^x) 


rij dC 


(6.19) 


where are components of the strain tensor. 

By expanding the expression, one can show that J(F,F) in Equation (6.19) is 
equivalent to 


J(F,F) = J(F) + J(F) + J X 


( 6 . 20 ) 


in which J(F) is the well known J-integral for stress intensity factor computation, 
J(F) is the ./-integral of the point force solutions, and 


Jr 



( Fik^ik T & ik&ik)d X j 




rij dC 


( 6 . 21 ) 


is the integral associated with the “cross-terms”. From the derivations in [16] and 
[67], we have: 


■J(F) = 0 and J x = 


77 

E 


( 6 . 22 ) 
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By rearranging Equation (6.20), we obtain a conservative integral for T-stress: 


T 


E 

1 


[J(F, F ) - J(F)] 


(6.23) 


In [67], the point force solution was obtained from an additional finite element 
calculation. This implies that for a specific problem, two finite element compu- 
tations need to be performed to obtain T-stress. In this study, we will use the 
analytical fields directly from the point force solution to evaluate T-stress. The 
displacement derivatives with respect to x and y for the point force auxiliary fields 
are: 


fJ'XjX — 
Uy,x = 


/ k cos 9 + cos 3 9 

4 n r 2 G 
f — ft sin 9 + sin 30 

4 7 rr 2 G 


(6.24) 

(6.25) 


U x ,y 

u y<y 


f 

— — — (1 + cos 2 9 + v cos” 9) sin 9 
Enr 
f 

— — (— sin 2 9 + v cos 2 9) cos 9 
Enr 


(6.26) 

(6.27) 


and the T-stress can be readily evaluated from Equation (6.23) using a finite 
element solution and analytical stress and displacement derivative fields of the 
point force. 


6.2.3 Analytical Equivalence between Betti’s Reciprocal 
and Eshelby Integrals 

To prove that both types of contour integrals for T-stress are analytically equiva- 
lent, we first observe that: 


^i,x 

a ij,x 


* 

U* 


a. 


i-3 


by setting the point force equal to one. Substituting this relationship into Equa- 
tion (6.14), we have: 


T 

E 


(<j*jiii — OijU*)rijdC 


c 


'fj *x^'i)r j d ( ' dC 


(6.28) 


From chain rules, we have (<JijU,) )X = (oy^rq) + and Equation (6.28) 

becomes: 


T 

E 


{p'ij'U'i) jx'ft'j dC dC dC 

C Jc 


(6.29) 



113 


By applying divergence theorem and recalling a lh j = 0 and e t j = - j + ) , 

the first term in Equation (6.29) can be expressed as: 


c 


ip ijV"i) ,x^j dC 


jj (Pij‘U , i),xj dA 


i&ik&ik^xj ) ^ j dC 


C 


c 


2 i&ik&ik T ^ikPik) d xj dC 


(6.30) 


Substituting Equation (6.30) into Equation (6.29), we then prove the analytical 
equivalence between Betti-Rayleigh reciprocal and Eshelby type of contour inte- 
grals. 


6.3 T-Stress Evaluation Using Finite Element 
Analyses 


6.3.1 Equivalent Domain Integral 

The conservative line integrals for T-stress derived above may not always be suit- 
able for use directly with results obtained by standard finite element methods. A 
procedure that converts a line integral into an equivalent area (or domain) integral 
by means of Gauss integral theorem is usually employed. The equivalent domain 
integral is known to have higher accuracy in extracting the desired integral values 
from given standard finite element solutions [78, 7]. 

Following Li et al. [78], for cracks in homogeneous bodies, the integrals in 
Equations (6.15) and (6.23) can be converted into their equivalent domain integrals 
over a closed area A. For the Betti-Rayleigh reciprocal type of contour integral, 
the counterpart is given by: 


T 



°ij u i E ) ( U dA 


For an Eshelby type contour integral, the domain integral is given by: 


(6.31) 


T 


J [J(F, F ) - J(F)} 

7 [If A i : ,1 7 + ,1 ’' + V, ‘ - ~ ' ^ FE + ^ <l : jdA 

~ If A H E uff - W FE S„] qjdAj 


(6.32) 



Figure 6.3: An equivalent domain integral. 


where the area A is a region as illustrated in Figure 6.3 and the function q is taken 
to be unity on , o and zero on , i. To quantify our computed results, we denote 
ro and ri as the distances from the crack tip to , o and , i along the 0 = 0 ray, 
respectively. The distance r x will be used to characterize the integration domain. 

6.3.2 Hierarchical p-version Finite Element Method 

To evaluate Equations (6.31) and (6.32) we compute displacements, strains and 
stresses using a hierarchical p-version finite element method. A comprehensive 
description of the discretization properties as well as implementation details of the 
“p-version” can be found e.g. in [141], We like to recall as an important property 
that, for problems with singularities, p-extensions converge with exponential rate 
if the mesh is properly refined towards the singularities. 

The element formulation used to obtain the results presented in Section 6.4 
restricts us to a standard, hierarchical polynomial basis for the finite element test 
and trial spaces. Although it is possible to introduce singular terms by the quarter- 
point mapping technique into hierarchical p-version formulations [109] these terms 
are less relevant for T-stress extraction. 

The use of high order p-version elements allows domain integrals to be com- 
puted during the postprocessing step on an integration mesh that can be completely 
independent from the finite element mesh. Examples for two different possibilities 
of postprocessing meshes are shown in Figure 6.4(a) and (b). If the solution is of 
sufficiently high quality so that jumps in the stresses are small, it is even possible to 
have elements of the integration mesh reaching over more than one element of the 
discretization without significant loss in accuracy of the domain integral. In this 
case, an additional step of locating the gauss points of the integration element in 
the corresponding finite element would become necessary but the implementation 
can be integrated much easier in a CAD environment. 

Numerical experiments show that it is not possible to reliably compute T-stress 
inside the elements directly adjacent to the crack tip. This is because, no matter 
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(a) independent (b) direct 

Figure 6.4: Postprocessing meshes for domain integral. 


how strong the mesh refinement towards the singularity is, the solution in the 
elements directly adjacent to the singularity can oscillate [111]. The oscillation 
behavior is restricted to the crack tip elements and to a lesser degree for the first 
layer of elements outside the crack tip elements. If the discretization error in the 
remaining domain is reduced sufficiently, then path independence can be observed 
in the numerical results. Therefore the domain integral can always be computed 
in the first or second layer outside the crack tip elements. 

To remove dependency of the obtained results from the error introduced by 
numerical integration we choose the number of integration points in each direction 
to be 15 for all values of p-extensions in all our computations. The additional 
computational cost for postprocessing of fracture parameters is orders of magni- 
tude less than the cost of the solution process and can be neglected. With this 
high integration order, numerical equivalence of the Betti-Rayleigh reciprocal and 
Eshelby type contour integrals can be directly observed. It may finally be noted 
that, using high order p-elements, it is very well possible to obtain high quality 
results for T-stress directly from contour integrals removing the need for this kind 
of postprocessing completely. 

6.3.3 Error Analysis and Accuracy Assessment 

In theory, the contour integrals developed herein are path independent. The finite 
element approximation, however, inevitably introduces discretization error. The 
quality of the obtained ./-integral or T-stress results does therefore depend on 
the location where the path independent integral has been evaluated. Thus, it is 
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important to quantify the error in computing T-stress using the path independent 
integrals and furthermore assess the accuracy of the obtained numerical results. 

For a given finite element model, the difference between the exact and finite 
element solutions is the discretization error. Thus, taking the Betti-Rayleigh re- 
ciprocal type contour integral for example, the error for T-stress, e T = \T — T FE \ 
at a certain integration path C is given as: 


e T = E\sf {a*jUi - (JijU* )rijdC— sf (a* jU EE - a EE u*)n j d,C \ 

Jc Jc 

Because of the different convergence rates for singular and non-singular terms 
in the solution, we shall observe that the discretization error is dominated mainly 
by the singular terms, if we extract our results in the singular-dominant zone. 
Thus we may postulate that the error of computed T-stress for Mode I problems 
(the effect of Mode II will be discussed later) in an integration domain rq away 
from the crack tip is: 


T - Kj 
e ~ e • — = 

Vo 

where e is a constant term related to the discretization error for a given discretiza- 
tion with a fixed polynomial degree of the shape functions. That is, the ratio of 
Ki/y/n can be factored out in e T similar to the asymptotic stress expansion. The 
relative error in T-stress 


(6.33) 


C-rei 


| T - T 


FE | 


K, 


T 


~ a 




(6.34) 


scales with Kj/(T y/r\). From Equation (6.34), we shall anticipate that the ac- 
curacy of the computed T-stress can be improved predominately by reducing the 
discretization error e, or by increasing the size of the integration zone (which may 
not always be practical). The assertion that the relative error in T-stress scales 
with the dimensionless parameter Kj /(Ty/r[) is also supported by the following 
observation: geometrically similar finite element models which differ only in scale 
(which implies that the integration path is likewise scaled) should give numerically 
identical error fractions in the computed T-stress (or any local stress measure- 
ment). 

We shall finally note that Equations (6.33) and (6.34) provide a powerful and 
useful measurement to study the accuracy of T-stress computations for finite el- 
ement analyses, as will be revealed in the following numerical results. One can 
certainly use the property of this error estimator to calibrate and regularize the 
computed T-stress for other standard, low-order p finite element codes, but this is 
beyond our discussion in this study. 
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(a) 


(b) 


Figure 6.5: Benchmark example for T-stress computation: (a) numerical model 
and (b) mesh with 6 layers of refinement (only 2 visible as shown) 


6.4 Numerical Results 

6.4.1 A Simple Benchmark Example 

In order to evaluate the accuracy of the proposed algorithm and its numerical 
outcome, it is desirable to compare it with an exact solution. Few exact solutions 
are known for the T-stress. For convenience, a problem with a simple geometry 
was chosen, with boundary conditions applied simulating arbitrary values of A/, 
A'//, and T. Since exact solutions are known, the problem may serve as a simple 
benchmark to study the accuracy of T-stress computations. 


6. 4. 1.1 Numerical Model 

As shown in Figure 6.5(a) an edge cracked square plate with a /L = 0.5 was mod- 
eled. Stresses obtained from the A'/, K n , and T related terms of the asymptotic 
expansion according to Equations (5.4), (5.5), and (5.6) were imposed as bound- 
ary conditions on the boundaries B 2 , B 3 , and B 4 , while the crack edges B x remain 
traction free. Thus the model represents a variant of an internal crack in an infinite 
plate under remote loading, such that A/, K n , and T can be varied independently. 
Because no body forces are present, the system is in equilibrium and the solution 
does exist. In addition to ensure uniqueness of the solution, displacement boundary 
conditions to prevent rigid body motions have to be prescribed (Figure 6.5(a)). A 
finite element mesh was constructed with six layers of refinement towards the crack 
tip resulting in a total of 52 elements (Figure 6.5(b)). The mesh was geometrically 
refined towards the crack tip with a progression factor of 0.15 as recommended in 
[MO]. 





6.4. 1.2 Numerical Results and Discussion 

For three different load combinations K t /T = 0.1, K r /T = 1.0 and K t /T = 10.0 
uniform p-extensions were performed, i.e. the polynomial degree p was changed for 
all elements uniformly between 1 and 11. The largest discretization with p = 11 
had 6169 degrees of freedom. In the postprocessing step of all computations the 
discretization error in energy norm ||e|[, the error in Kj and the error in T were 
computed. 

The first observation to be made from the results is that the T-stress values 
computed from Betti-Rayleigh reciprocal and Eshelby type integrals coincide for 
all computations and all degrees of p up to machine accuracy. The equivalence 
proven in Section 6.2.3 is therefore also visible in the numerical results, if the 
integration order for the domain integral is sufficiently high. 

In Table 6.1 values for the J-domain integral, A'/, and T-stress, computed on 
four different integration domains are tabulated. The domains are coincident with 
the first to the fourth layer of elements surrounding the crack tip elements. The 
load parameter Kj/T was 1.0 and the polynomial degree of the elements was p = 6 
corresponding to 1899 degrees of freedom. Since the discretization error is virtually 
non-existent for this case, path independence of Kj obtained from J as well as the 
T-stress term can be observed. 

In Figure 6.6(a)-(c) convergence curves of computed T-stress values are shown 
for various Kj/T combinations. In each figure the error in T-stress computed on 
the first layer (LI, rq = 7.59375 x 10 -5 ) and the second layer ( L2 , r x = 5.0625 x 
10 -4 ) of elements outside the crack tip is plotted over the number of degrees of 
freedom. Each mark indicates a polynomial degree p. As an indication of the global 
convergence behavior of the solution, the error in energy norm is also plotted for 
all load combinations. While the energy norm curves show the typical S-shape (i.e. 
exponential convergence rates until p equals the number of refinement layers), the 
curves for the T-stress show exponential convergence rates throughout the entire 
p-range. It is further observed that, for constant rq, the relative error in T-stress 
increases with Ki/T but, as apparent in Table 6.2, the convergence rate is exactly 
the same for all loading values as we shall expect from Equation (6.34). Finally 
it is also observed that curves for computation of T-stress on the second layer 
outside the crack tip elements are slightly smoother, especially for low orders of 
p, indicating that oscillatory behavior of the solution is restricted to the crack tip 
and first layer of elements. 

Table 6.2 summarizes the computed errors of A7, T and the energy norm on 
the first layer of elements away from the crack tip for all three Kj/T ratios. From 
these values it becomes apparent that the relative error in K t is independent of 
Kj/T. The convergence rate of the relative error in energy norm calculated at 
p = 4 and p = 8 from 
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Table 6.1: Numerical Path Independence of J, Aj, and T for p = 6 and Kj/T = 1 
for Benchmark Example 


Path 

r i 

J-domain 

K, 

T 

4 

2.25 x 10- 2 

0.91012 

1.0001 

0.9997 

3 

3.375 x 10- 3 

0.91019 

1.0001 

0.9993 

2 

5.0625 x 10“ 4 

0.91026 

1.0001 

0.9980 

1 

7.59375 x 10“ 5 

0.91032 

1.0002 

0.9942 


*plane strain problems with E = 1.0 and u = 0.3 


is (3 = —2.15 for KijT = 1 which is about the same as /3 = —2.2 for K/T = 10. 

We further evaluated the model with various values of A'j, T, and rq for var- 
ious p-extensions to obtain the results plotted in Figure 6.7. Both the order of 
p and Ki/(Tyf¥i) are seen to have a significant effect on the accuracy of the so- 
lution. The results are expected from Equation (6.34), since the relative error in 
T-stress is proportional to the discretization error e and the dimensionless parame- 
ter KiKT-yJrl). We finally note that with p = 11, a relative error in the calculated 
T-stress approaching 1 x 10 -6 % is achieved. 

In the process of simulations, various Kjj/Kj ratios were also evaluated, and 
the effect of K n was found to be negligible for values up to K n /Kj = 100. Thus 
the effect of K u was not given further consideration in this study. 

From the above we observe that for problems where T-stress is small com- 
pared to A'/, it is more difficult to accurately evaluate. Perhaps most significantly, 
the above exercise identifies the conditions under which we may with confidence 
calculate T-stress with very high accuracy. 


6.4.2 Fracture Specimens 

In this section, the numerical results for various fracture specimens are evaluated. 
This serves two purposes: one is to demonstrate that the convergence and accu- 
racy of computed T-stress can be observed easily with the p-extensions, and the 
other is to compare our results for various fracture specimen configurations with 
numerical values from the literature. We detail our results for double cantilever 
beam (DCB) specimens due to its practical importance in obtaining an accurate 
T-stress to characterize crack turning behavior [103]. For other fracture specimens, 
we tabulate our computed T-stress for comparison. 

6. 4. 2.1 Double Cantilever Beam (DCB) Specimen 

DCB fracture specimens are known to have large positive T-stress that may cause 
crack path instability under pure Mode I conditions [46, 26]. In order to compare 
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Table 6.2: Relative Error of Kj and T Computed at Integration Domain on First 
Layer Away From Crack Tip and Relative Error in Energy Norm for 
Benchmark Example 


n = 7.59375 x 10“ 5 , K t /T = 1/10 

P 

DOFs 

relative error (%) 

K, 

T 

INI 

4 

919 

0.456 

0.7693 

0.0781 

8 

3295 

4.3810 x 10“ 4 

2.7548 x 10~ 3 

0.004 

11 

6169 

8.1660 x 10~ 5 

3.0 x 10" 5 

0.002 


ry = 7.59375 x 10- 5 , K,jT = 1/1 

P 

DOFs 

relative error (%) 

K, 

T 

e|| 

4 

919 

0.456 

7.693 

0.656 

8 

3295 

4.3810 x 10" 4 

2.7548 x 10“ 2 

0.042 

11 

6169 

8.1650 x 10~ 5 

2.9995 x 10" 4 

0.022 


ri = 7.59375 x IQ' 5 , K r /T = 10/1 

P 

DOFs 

relative error (%) 

K, 

T 

|e|| 

4 

919 

0.456 

76.93 

1.4675 

8 

3295 

4.3810 x 10" 4 

0.27548 

0.073 

11 

6169 

8.1660 x 10" 5 

2.996 x 10- 3 

0.046 


^results are independent of the length units associated with r± and Kj/T 





relative error, e^ pJ x 100 j(%|) 
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p=2 3 4 5 6 7 8 91011 p=3 4 5 6 7 8 910 11 



(a) K,jT = 1/10 (b) Ki/T = 1/1 



(c) K,/T = 10/1 


Figure 6.6: Convergence of T-stress and energy norm for the benchmark example. 

T-stress values are computed on the first layer (LI, r± = 7.59375 x 10 -5 ) 
and the second layer (L2, r\ = 5.0625 x 10 -4 ) of elements outside the 
crack tip. 
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(a) (b) 

Figure 6.7: Accuracy assessment with p-extension for the benchmark example: (a) 
the relative error of calculated T-stress, ej d , versus dimensionless pa- 
rameter Ki/(Ty/ri), and (b) the ratio of ej d over A'//(T v /r7) versus 
the polynomial degree of the shape functions, p. 

numerically computed T-stress values to those found in the literature, a DCB 
configuration shown in Figure 6.8(a) with h/w = 0.2 and a/w = 0.5 was modeled. 
Again a finite element mesh with six layers of refinement towards the crack tip was 
constructed, as shown in Figure 6.8(b). To eliminate any influence of perturbations 
from point forces the load was introduced as distributed forces along one half of 
the hole edges. 

With the mesh fixed, a p-extension was performed, i.e. p was increased between 
1 < p < 11. The discretization at p = 11 had 9277 degrees of freedom and, because 
no analytical solutions are available for this specimen, it was used as a reference 
solution. K] and T were computed in two different domains corresponding to the 
first and second layer of elements away from the crack tip with ri/a = 2.025 x 10~ 4 
and ri/a = 1.35 x 10 -3 . We note that the value of the error parameter K I /{T- s Jr[) 
for the integration domain at the first layer is 42.2; thus the relative error of the 
reference solution is estimated to be on the order of 10 -5 % based on Figure 6.7. 

The convergence of the computed normalized I\ I and T-stress values can be 
observed in Figure 6.9 and Figure 6.10, respectively, where for each polynomial 
degree from p = 3 to p = 8, A'/ and T-stress computed on the two integration 
domains are plotted. For engineering accuracy, i.e. to obtain a relative error of 1 
percent in T-stress, a degree of p = 6 corresponding to 2917 degrees of freedom 
is necessary. Again path independence of the computed integrals can be observed 
from p = 5 and up. 

To compare our results with values from different sources published in the 
literature, a normalized stress biaxial parameter B defined in [77] is introduced 
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Table 6.3: Computed Values of K h T, and B for the DCB Specimen 


DCB (a/w = 

Od" 

o 

II 

p-a 

o 

II 


SOURCES 

Ki 

(jy/i ra 

T 

a 

D TV 7ra 

B = K, 

Present (p = 11) 

3.9225 

11.5745 

2.9508 

Leevers and Radon [77] 

- 

- 

2.942 

Cardew et al. [16] 

- 

- 

2.829 

Ivfouri [67] 

- 

- 

2.956 

Fett [45] 

3.9307 

11.5304 

2.933 


where 


_ T s/tox 
= Ki 


(6.35) 


Table 6.3 summarizes the normalized stress intensity factor Kj/ ra), the nor- 
malized T-stress, and B with p = 11 as well as computed values from [77, 16, 67, 
45]. A difference of up to 4.2% in the results for B is observed in comparison with 
the sources from the literature. 


6. 4. 2. 2 T-stress for Various Fracture Specimen Configurations 

Computed T-stress in various fracture specimen configurations including middle 
crack tension specimen (MT) and single edge notch specimen subjected to tension 
(SENT) and pure bending (SENB) as shown in Figure 6.11 was evaluated. All the 
meshes were constructed with six layers of refinement towards the crack tip and a 
progression factor of 0.15. 

Table 6.4 summarizes our computed results for K r , T- stress, and B. The results 
are compared with values from different sources published in the literature. 


6.4.3 Discussion: Numerical Assessment of T-stress Com- 
putation Using a p-version Finite Element Method 

Throughout the numerical examples, we have demonstrated that using the path 
independent integrals with hierarchical, p- and hp-version finite element methods 
proves to be a powerful tool to obtain highly accurate numerical results for T- 
stress. The convergence and accuracy of computed T-stress values are observed 
easily and confidently with the p-extensions in the benchmark example of known 
exact solutions, and the error correlates reliably with the error estimator, Kj/ \fr{. 
It is thus inferred that the results presented for the fracture specimen geometries 
are of comparable accuracy to the benchmark at equal values of K / / ^/rf, and are 
thus numerically exact to the significant digits given in the tables. 




(b) 


Figure 6.8: (a) Double cantilever beam (DCB) specimen configuration, and (b) a 
hp-version finite element model for DCB specimen with 6 layer refine- 
ment (only 2 visible). 
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p = 3 4 5 6 7 8 



Figure 6.9: Convergence of normalized Kj for p-extension for the DCB specimen. 


p = 3 4 5 6 7 8 



Figure 6.10: Convergence of normalized T-stress for p-extension for the DCB spec- 


imen. 
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Table 6.4: Computed Values of Ki, T, and B for MT, SENT, and SENB Specimens 


SOURCES 

K, 

(Jx/jra 

T 

a 

75 Tx/na, 

B = K, 

MT (2a/w = 0.3, h/w = 1.0) 

Present (p = 11) 

1.1232 

-1.15536 

-1.0286 

Leevers and Radon [77] 

- 

- 

-1.0255 

Cardew et al. [16] 

- 

- 


Fett [45] 

- 

-1.1557 

-1.028 

Isida [61] 

1.123 

- 

- 

SENT (a/w = 0.3, h/w = 12) 

Present (p = 11) 

1.6598 

-0.61033 

-0.36771 

Sham [123] 

1.6570 

-0.61425 

-0.37070 

SENT (a/w = 0.5, h/w = 12) 

Present (p = 11) 

2.8246 

-0.42168 

-0.14929 

Sham [123] 

2.8210 

-0.43142 

-0.15293 

SENB (a/w = 0.3, h/w = 12) 

Present (p = 11) 

1.1241 

-0.079177 

-0.070436 

Sham [123] 

1.1220 

-0.082404 

-0.073444 

SENB (a/w = 0.5, h/w = 12) 

Present (p = 11) 

1.4972 

0.39749 

0.26549 

Sham [123] 

1.4951 

0.39112 

0.26160 
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a M <r = 6M/t 2 

(a) MT (b) SENT (c) SENB 


Figure 6.11: Various fracture specimen configurations used for T-stress computa- 
tion. 

6.5 FRANC3D/STAGS Results 

Both the Betti-Rayleigh reciprocal and Eshelby types of domain integrals were 
also implemented in the FRANC3D/STAGS software program. We note that the 
polynomial degree of shape functions for the quadrilateral shell element used in 
FRANC3D/STAGS is the lowest, i.e., p = 1. Thus, to obtain an acceptable 
accuracy for T-stress, we need to introduce a highly focused mesh near the crack 
tip and/or evaluate the domain integral sufficiently away from the crack tip. 

6.5.1 Two-Dimensional Problems 

The DCB specimen was studied using FRANC3D/STAGS. An all-quadrilateral 
element meshing algorithm developed by Potyondy et al. [107] was used to generate 
a graded mesh with a high mesh density near the crack tip and coarser away from 
the crack tip. Figure 6.12 shows a graded finite element mesh with a crack tip 
template. Computed values of Kj and T evaluated at the third layer away from 
the crack tip (rq/a = 0.02) are within 0.8 and 1.2% of the reference solutions in 
Table 6.3, respectively. 
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Figure 6.12: A focused finite element mesh used in FRANC3D/STAGS. 

6.5.2 Thin Shell Problems 

The methods developed herein are mainly for two-dimensional linear elastic prob- 
lems. Further study is needed to derive its Lagrangian counterpart for shell struc- 
tures subjected to large displacements and rotations. 

6.6 Summary 

Two types of path independent integrals for T-stress computations, one based 
on the Betti-Rayleigh reciprocal theorem and the other based on Eshelby’s en- 
ergy momentum tensor were studied. Analytical as well as numerical equivalence 
between the two integrals was found. To quantify and assess the accuracy of com- 
puted values, a novel error estimator for T-stress was proposed. Specifically, it 
was found that the error of the computed T-stress is proportional to the ratio of 
stress intensity factor divided by the square root of the characteristic dimension of 
the integration domain where the path independent integral is evaluated. Using a 
highly accurate hierarchical p- and hp-version finite element code, the convergence 
and accuracy of computed values were observed easily and confidently, and the er- 
ror of the computed T-stress correlated reliably with the proposed error estimator. 
We conclude that the path independent integrals, in conjunction with hierarchical, 
p- and hp-version finite element methods, provide a powerful tool to obtain highly 
accurate numerical results for T-stress. 

We then evaluated numerical results using the FRANC3D/STAGS program. 
Because the polynomial degree of the shape functions of the shell elements was the 
lowest (i.e., p = 1), a highly focused mesh near the crack tip and a remote inte- 
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gration path were needed to extract the T-stress from the finite element solutions 
without losing too much accuracy. Numerical results showed that the errors in T- 
stress for DCB specimens could be kept well below 2% using FRANC3D/STAGS. 



Chapter 7 

Curvilinear Crack Growth 
Simulations 


Curvilinear crack growth in coupon tests and in full-scale curved panel tests are 
analyzed in this chapter. Crack growth direction is predicted using the directional 
criteria developed in Chapter 5. The predicted crack trajectories are compared 
with those observed in the tests. 


7.1 Curvilinear Crack Growth Simulations For 
DCB Fracture Specimens 

Simulations of crack growth in double cantilever beam (DCB) fracture specimens 
were performed. The predicted crack trajectories were compared with the exper- 
imental measurements. Among all the possible parameters that could affect the 
crack trajectory prediction, only T-stress, r c , fracture toughness orthotropy, and 
the length of the crack growth increment were examined. 

7.1.1 Description of Experiment 

DCB specimens made of 0.09 inch thick, 2024-T3 aluminum alloy were tested at 
the NASA Langley Research Center in cooperation with the McDonnell Douglas 
Company (now Boeing). A brief description of the tests is presented below. More 
information about the fracture tests can be found in [103, 101, 104], 

The dimensions and material properties of the test specimens are shown in Fig- 
ure 7.1. Stable crack growth in the L-T and T-L orientation under a monotonically 
increasing load was conducted. Fatigue crack growth in the L-T orientation under 
a low stress level of cyclic loading was performed. The test matrix is summarized 
in Table 7.1. Only one test per configuration was performed. 

The final cracked configurations in the L-T orientation are shown in Figure 7.2. 
The crack path observed in stable crack growth was different from that in fatigue 
crack propagation. For specimens under stable tearing, the crack turned sharply 
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Figure 7.1: The DCB specimen configuration. 


Table 7.1: Test Matrix for DCB Specimens 


Specimen ID 

Initial Crack (in.) 

Type of Crack Growth 

2024LT-4 

5.507 

stable tearing 

2024TL-5 

5.47 

stable tearing 

2024LT-6 

5.104 

fatigue crack growth 




( a ) 


(b) 


Figure 7.2: (a) Stable tearing (2024LT-4) and (b) fatigue crack growth (2024LT-6) 
in the L-T orientation observed in DCB specimens (after [101]). 

away from its initial crack tip. For specimens subjected to fatigue loading, a much 
smoother crack path was observed. 

7.1.2 Description of Simulation 

Curvilinear crack growth analyses were conducted using the FRANC3D/STAGS 
software program. As described in Section 1.3, to simulate crack growth where 
the crack trajectory is not known a priori, continual updating of the geometry is 
required. The procedure of simulating crack growth consists of the following steps 
[106, 108]: 

step 1: generate a STAGS finite element model from FRANC3D. 

step 2: obtain the equilibrium state by executing the STAGS code. 

step 3: compute the fracture parameters and determine the direction of crack 
growth. 

step 4: decide the amount of crack growth extension and propagate the crack. 

The process is repeated until a suitable termination condition is reached. The 
crack growth alters the geometric model in FRANC3D and leads to localized mesh 
deletion. The deletion region is remeshed automatically using an all-quadrilateral 
element meshing algorithm [107]. 
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Stress intensity factors, K t and K n , were computed using the equivalent do- 
main J integral with the mode-separation method [15]. T-stress was computed 
using the equivalent domain integral method derived from the Betti-Rayleigh re- 
ciprocal theorem (Equation (6.15)). Several integration paths about the crack tip 
were evaluated to assure the accuracy of the computations. The propagation angle 
was predicted based on the maximum tangential stress theory developed in Chap- 
ter 5, that is, Equations (5.17) and (5.24) for the isotropic and orthotropic media, 
respectively. We assumed that both fatigue and stable crack growth in the DCB 
specimens can be analyzed using linear elastic fracture mechanics (LEFM). This 
assumption may not hold for stable crack growth. We will further discuss this 
issue in Section 7.3. 

7.1.3 Numerical Results 
7. 1.3.1 Effect of T-stress and r c 

The specimens under stable crack growth were analyzed first. The crack was grown 
in 0.2 inch increments. Crack growth direction was predicted by the isotropic di- 
rectional criterion, i.e., Equation (5.17). Figure 7.3 depicts four predicted crack 
trajectories with various magnitudes of r c and the experimental measurements. 
Figure 7.4 shows the computed deformed shapes and the corresponding finite ele- 
ment meshes used in curvilinear crack growth simulations. Predicted and measured 
results shown in Figure 7.3 indicate that: 

1 . The predicted crack path for r c = 0 coincides with the straight line ahead of 
the initial crack. For this special case, we note that the directional criterion 
reduces to the Erdogan and Sih’s criterion [41]. 

2. For r c = 0.05 inch, the predicted crack path initially follows a zigzag along 
the straight path, but deviates from it at about 1.6 inches of crack extension. 

3. For r c > 0.06 inch, predicted crack paths turn sharply away from the initial 
crack tip. 

4. The case with r c = 0.09 inch best correlates the experimental data and 
numerical results for stable tearing. 

5. Notable difference between the measured crack paths in the L-T and T-L 
orientations are observed. 

Similar trends were observed for fatigue crack growth. Figure 7.5 plots three 
predicted crack trajectories with various magnitudes of r c and the experimental 
measurements. The case with r c = 0.06 inch best correlates the experimental data 
and numerical results for fatigue crack growth. 
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o 2024LT-3 (side 1) 


• 2024LT-3 (side 2) 

□ 2024TL-4 (side 1) 




(b) 

Figure 7.3: Predicted and measured crack trajectories for DCB specimen under 
stable crack growth: (a) overall crack trajectories, and (b) crack paths 
in the focused region (Equation (5.17) with various magnitudes of r c ; 
A a = 0.2 in.). 
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(a) A a = 0.6 inch 



(b) A a = 2.0 inch 

Figure 7.4: Computed deformed shapes and the corresponding finite element 
meshes used in the curvilinear crack growth simulations in DCB spec- 
imens (isotropic case with r c = 0.09 inch). 
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o 2024LT-6 (side 1) 
• 2024LT-6 (side 2) 

(in.) 



Figure 7.5: Predicted and measured crack trajectories for DCB specimen under 
fatigue crack growth (Equation (5.17) with various magnitudes of r c ; 
A a = 0.2 in.). 

7. 1.3. 2 Length of Crack Growth Increment 

For LEFM problems, no adaptive scheme is currently implemented to control the 
length of curvilinear crack growth increment. Instead, during the simulation pro- 
cess, the analyst needs to decide and specify the amount of crack extension at 
the crack tip. To investigate possible effects of the crack growth increment on 
crack trajectory prediction, the crack growth simulation was performed again with 
A a = 0.1 inch for r c = 0.09 inch. The length of crack growth increment has a 
minor effect on the crack trajectory prediction as shown in Figure 7.6. 

7. 1.3. 3 Fracture Toughness Orthotropy 

From the laboratory observation, the specimen in the L-T orientation turned 
sharper than that in the T-L orientation under stable tearing. This is thought 
to be related to the possible difference between the fracture resistance along the 
transverse (T) direction compared to that along the longitudinal (L) direction. 
In subsequent analyses, a simple elliptical function presented in Section 5.3.4 was 
used to incorporate the effect of fracture orthotropy. The fracture toughness was 
assumed to be 10% higher in the T than in the L direction. As shown in Fig- 
ure 7.7, the predicted crack trajectory in the L-T orientation agrees better with 
the experimental measurements than the isotropic prediction. However, the pre- 
dicted trajectory in the T-L orientation deviates from the observed crack path. 
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o 2024LT-3 (side 1) 

• 2024LT-3 (side 2) 

□ 2024TL-4 (side 1) 

■ 2024TL-4 (side 2) 



Figure 7.6: Predicted and measured crack trajectories for DCB specimen under 
stable crack growth: effect of the length of crack growth increment 
(Equation (5.17) with r c = 0.09 inch). 

7.1.4 Concluding Remarks for Curvilinear Crack Growth 
Simulation in DCB Specimen 

Curvilinear crack growth in thin, metallic DCB specimens was studied. For this 
specific configuration, cracks showed a tendency to turn away from the initial 
crack tip under pure Mode 1 conditions. The crack growth directional criterion, 
incorporating the T-stress effect, was capable of capturing the essence of crack 
turning under such circumstance. The predicted results were in good agreement 
with the experimental measurements. 

The fracture toughness orthotropy was introduced to explain the difference be- 
tween the measured crack paths in the L-T and T-L orientations. The orthotropic 
directional criterion showed some promise to correlate the experimental data, but 
some disagreement between predicted and measured results was observed for the 
specimen in the T-L orientation. One possible explanation is that the magnitude of 
r c along the T direction may be different than that along the L direction. Incorpo- 
rating different magnitudes of r c along T and L directions could certainly provide 
better correlation with experimental results. But additional tests and analyses 
need to be conducted to justify this assertion. 



138 


o 2024LT-3 (side 1) 

• 2024LT-3 (side 2) 

□ 2024TL-4 (side 1) 

■ 2024TL-4 (side 2) 



(in.) 


Figure 7.7: Predicted and measured crack trajectories for DCB specimen under 
stable crack growth: effect of the fracture toughness orthotropy (Equa- 
tion (5.24) with r c = 0.09 inch; = 1.1; A a = 0.2 in.). 




139 


Material: Skin and tear strap, 2024 13 clad, \ . 10,500 ksi, v 033 

stringer and frame, 7075 TO clad. It 10,700 ksi, v 0.33; 


Radius: 74 inches 
Skin: 0.036 in. thick 
Tear strap: 0.036 in. thick 
Stringer: 0.028-in. -thick hat section 
9.25-in. spacing 

Frame: 0.040-in. -thick Z-section 
20-in spacing 


Tear straps 



Stringer 


Stringer clip* 


Skin 


Frame 


Tear strap 


Frame 


Figure 7.8: Structural features of a narrow body fuselage panel (modified after 
[85]). 

7.2 Curvilinear Crack Growth Simulations For 
Fuselage Structures 

Simulations of curvilinear crack growth in a generic narrow body fuselage panel 
were performed. The predicted crack trajectories were compared with the mea- 
sured values from a full-scale pressurization test. The problem demonstrates the 
applicability of the direction criteria developed herein for predicting curvilinear 
crack growth in fuselage structures. 

7.2.1 Description of Experiment 

A narrow body fuselage panel with tear straps, stringers, stringer clips, and frames 
was tested by the Boeing Commercial Airplane Group. Skins and tear straps were 
made of 0.036 inch thick, 2024-T3 clad aluminum alloy. Stringers, frames, and 
stringer clips were made of 7075-T6 clad aluminum alloy. The tear straps were 
hot-bonded to the skins at midbay and at each frame station. The structural 
features of the test panel are shown in Figure 7.8. More information about panel 
dimensions can be found in [106, 48]. 

The panel had a 5.0 inch initial saw cut in the T-L orientation centered on 
the midbay tear strap and just above the stringer tear strap. The saw cut went 
completely through both the skin and midbay tear strap. The panel was inserted 
into a test fixture with a radius of curvature of 74 inches to match narrow body 
airplanes. A cyclic pressure of 7.8 psi was applied to propagate the crack. During 
the test, the positions of the crack tips were recorded. The detailed test data can 
be found in [106]. 



Figure 7.9: Finite element model for the narrow body fuselage panel. 


7.2.2 Numerical Model 

The entire curvilinear crack growth simulation consists of more than 20 inches of 
crack extension. As a result, using a global-local hierarchical modeling approach 
could require continual updating of the boundary conditions from the preceding 
model in the hierarchy due to the crack growth. This would increase efforts sub- 
stantially in performing the numerical analyses. For this specific problem, only 
internal pressure was applied to the structure, thus a simple numerical model us- 
ing symmetric boundary conditions might suffice to simulate the panel test. 

In this study, a 4-stringer-bay wide and 2-frame-bay long panel was analyzed. 
All structural components including skins, stringers, and frames were modeled 
by quadrilateral shell elements. Each node of the shell element has six degrees 
of freedom. A typical finite element mesh used in the simulation is shown in 
Figure 7.9. 

Geometrically nonlinear analyses were performed. Pressure loading was applied 
on the skin of the shell model. Symmetric boundary conditions were imposed on 
all the boundary edges of the model to simulate a cylinder-like fuselage structure. 
Uniform axial expansion was allowed at one longitudinal end. On this boundary 
edge, an axial force equal to (PR/2) ■ L was assigned where P is the applied 
pressure, R is the radius of the panel, and L is the arc-length of the edge. 

7.2.3 Fracture Parameter Evaluation 

Deformation and stress fields near the crack tip were used to compute fracture pa- 
rameters for crack growth simulations. The modified crack closure integral method 
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was used to compute the membrane and bending stress intensity factors (K T , K n , 
ki, fc 2 ) [106, 142], The crack growth directional criteria, Equations (5.17) and 
(5.24) for the isotropic and orthotropic media, were used to predict the propaga- 
tion angle in thin shell structures. 

The equivalent domain integral method for T-stress developed in Chapter 6 
is only valid for two-dimensional problems. The derivation of its counterpart for 
shell structures subjected to large displacements and rotations is not yet available. 
Instead, a simple displacement correlation method was used to evaluate the T- 
stress [69, 134, 104], 

7.2.4 Numerical Results 

7.2.4. 1 Effect of T-stress and r c 

The effect of T-stress and r c on crack trajectory prediction was studied first. 
Crack growth direction was predicted by the isotropic directional criterion (Equa- 
tion (5.17)). Figure 7.10 plots the predicted crack trajectories with r c = 0 and 
r c = 0.09 inch as well as the experimental measurements. Figure 7.11 shows the 
computed deformed shapes during curvilinear crack growth. Bulging caused by 
the applied pressure is observed. Moreover, severe flapping is predicted as the 
crack turns. Figure 7.12 shows the computed stress intensity factors and T-stress 
versus the half crack extension at the right crack tip. The sign conventions of 
stress intensity factors follow those in [106]. Predicted results suggest: 

1. The T-stress has a very mild influence on the early crack trajectory prediction 
because of its small magnitude. But as the crack approaches the tear strap, 
T-stress increases and plays an important role in the crack turning prediction. 
For the case with r c = 0.09 inch, a sharp turning caused by T-stress is 
predicted as the crack approaches the tear strap. 

2. The computed fracture parameters for r c = 0 and r c = 0.09 inch are compa- 
rable at the early stage of curvilinear crack growth. However, sharp turning 
as the crack approaches the tear strap alters the deformation and stress fields. 
This drastically changes the computed values of fracture parameters. 

3. Predicted crack paths from both numerical simulations at the right and left 
crack tips are almost symmetric about the midbay, but the measured crack 
paths are not. This observation gives a preliminary indication of the experi- 
mental scatter that might occur in the panel test. 

7. 2. 4. 2 Effect of Fracture Toughness Orthotropy 

The predicted crack growth trajectories depicted in Figure 7.10 are comparable 
to the experimental measurements, but with some discrepancy. The disagreement 
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a — predicted (r c = 0) 

o— predicted (r c = 0.09 inch) 
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Figure 7.12: Computed stress intensity factors and T-stress versus half crack ex- 
tension. The hollow and solid markers denote the computed fracture 
parameters for the isotropic case with r c = 0 and r c = 0.09 inch, 
respectively. 
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frame line 


— measured 

-o- predicted (isotropic) _ 

—a— predicted (orthotropic, = 1-1) 

midbay frame line 



Figure 7.13: Comparisons between predicted and measured crack trajectories 
(isotropic and orthotropic cases with r c = 0.09 inch). 


during early stages of crack growth might be related to the fracture toughness 
orthotropy of the fuselage skins. 

In subsequent analyses, the orthotropic directional criterion, i.e., Equation (5.24), 
was used to predict the propagation angle. From the coupon test res ults, the frac- 
ture toughness for this material and thickness was about 100 ksivdnch in the L 
direction and 105-120 ksiydnch in the T direction [106]. Thus, the fracture tough- 
ness was assumed to be 10% higher in the T than in the L direction. The predicted 
crack trajectories with r c = 0.09 inch were compared with those from the isotropic 
prediction and experimental measurements. As shown in Figure 7.13, during early 
stages of crack growth, the predicted trajectories for the orthotropic case agree 
better with the experimental measurements than the isotropic case. Crack growth 
simulation with fracture orthotropy also predicts crack turning as the crack ap- 
proaches the tear strap. Yet, when the crack grows further into the tear strap 
region, the inclusion of fracture orthotropy adversely alters the crack path predic- 
tion and does not predict flapping as observed in the panel test. 

Several possible reasons may explain why the current methodology including 
the fracture toughness orthotropy does not predict the desired flapping and should 
be examined in future research: 

1. Characteristic feature of fracture orthotropy in the tear strap region — The 
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material orientation in the tear strap differs from that in the skin (?'. e., the 
transverse direction in the tear strap is along the longitudinal direction of 
the skin and vice versa). As a result, the material characteristics in this over- 
lapped region may behave like a quasi-isotropic material with less fracture 
toughness orthotropy. 

2. Occurrence of debonding between the skin and tear strap — In the current 
model, we assume that the skin and tear strap are perfectly bonded. However, 
as the crack grows into this region, the adhesive bond between the skin and 
tear strap is likely to fail. This inevitably alters the local crack tip stress 
fields and would consequently affect the crack growth behavior. 

3. Thin-shell representation of three-dimensional behavior — The thin-shell ap- 
proximation does not capture all the three-dimensional complexities of the 
problem in the vicinity of the tear strap, particularly in the crack tip region. 
Further study on three-dimensional crack growth simulations is needed to 
quantify the three-dimensional effect on crack turning prediction. 

4. Sources of error from the computed fracture parameters — Accurate stress in- 
tensity factor and T-stress evaluations in this region are crucial to predict 
crack turning. Current crack growth simulations use a low-order polyno- 
mial degree of shape functions for thin-shell finite element analyses and use 
a displacement correlation method to extract the T-stress term from the fi- 
nite element solutions. Further study using adaptive and higher order shell 
finite element analyses would improve the accuracy of numerical computa- 
tion. Other numerical methods, for example, path independent integrals for 
geometrically nonlinear shells would also improve the accuracy of fracture 
parameter evaluations. 

5. Validity of the LEFM approach — The crack growth directional criterion and 
its subsequent curvilinear crack growth simulations explicitly assume that the 
crack is grown under small scale yielding conditions. Yet, as the length of 
the fatigue crack extends to a sufficient amount, stable tearing and extensive 
plasticity are likely to occur. The active plastic zone and accumulated plastic 
wake due to stable tearing would likely affect the crack growth prediction. 

6. Validity of the magnitude of r c — In the current study, the parameter r c is as- 
sumed to be a constant magnitude during the entire curvilinear crack growth 
simulations. Also the same constant magnitude of r c is used in the T and L 
directions. The magnitude of r c used in the current simulation (0.09 inch) 
is mainly based on the predicted results that best correlate the crack tra- 
jectories observed in the DCB tests and the isotropic results that predict 
crack flapping in the current fuselage model. Further study on the appro- 
priate experimental methods and numerical simulations for determining the 
magnitude of r c is needed to validate the current approach. 
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7.2.5 Comparisons with Previous Studies 

Potyondy et al. [106, 108] and Chen et al. [21] have reported numerical simulations 
for this problem previously. Both studies analyzed early curvilinear crack growth 
but did not address the issue of sharp turning as the crack approaches the tear 
strap. It is, nevertheless, of interest to compare these results with the current 
prediction. This serves two purposes: one is to assess the accuracy of the simulation 
through comparisons with independent numerical simulations and the other is 
to show alternate modeling representations that may affect the crack trajectory 
prediction. 

In Potyondy et al. [106, 108] and Chen et al. [21], a global-local hierarchical 
modeling approach was used to model the panel test. Three hierarchical modeling 
levels were employed, comprised of a global shell model, a 6x6 bay stiffened panel 
model, and a 2x2 bay stiffened panel model. Crack growth was only performed 
in the 2x2 bay model, the lowest level in the hierarchy. The kinematic boundary 
conditions on the 2x2 bay model were not updated during crack growth. Also, the 
boundary conditions applied to the global shell model corresponded with an open 
cylinder. Thus, the longitudinal stress in this numerical model is expected to be 
less than that in the test fixture, since the test fixture is a closed cylinder. 

The directional criterion used in Potyondy et al. and Chen et al. corresponds 
to the Erdogan and Sih directional criterion [41], i.e., Equation (5.18); thus com- 
parisons are made with the isotropic prediction with r c = 0. Figure 7.14 shows the 
predicted crack growth trajectories from previous and current studies as well as 
experimental measurements. We note that the initial crack location in Potyondy’s 
simulations was modeled at 0.45 inch away from the intersection of the skin and 
stringer due to limitations in the previous version of the FRANC3D program. Fig- 
ures 7.15 and 7.16 show the computed stress intensity factors at the right crack 
tip in comparison with [21] and [106, 108], respectively. From these results, we 
conclude: 

• The applied axial force used to model the longitudinal stress caused by a 
closed cylinder has little influence on the computed stress intensity factors. 
This can be seen from the computed values shown in Figure 7.15 at zero 
crack extension; two numerical simulations at this stage basically represent 
the same boundary conditions and crack configuration except an axial force 
was applied in the current model. 

• The fact that the kinematic boundary conditions were not altered during 
crack growth in the previous studies has a mild affect on the crack trajectory 
prediction and stress intensity factor computation. In the previous studies, 
the kinematic boundary conditions used in the lowest level in the hierarchy 
were obtained from a global model with an initial 5.0 inch crack. We can 
certainly conclude that the driving force for this case would be less than the 
one with updated boundary conditions as the crack grows. This is properly 
reflected on the computed K r values shown both in Figures 7.15 and 7.16. 
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— measured 

-o- predicted (current, r c = 0) 

—a— predicted (Chen et al. [21]) 

— o— predicted (Potyondy et al. [106, 108]) 


frame line midbay frame line 



Figure 7.14: Comparisons between predicted and measured crack trajectories 
(isotropic directional criterion with various magnitudes of r c ). 

The issue seems to have little effect on the computed values of K n , since 
they remain more or less the same for all cases. This leads to a lower ratio of 
Kn/K / in the current model with updated boundary conditions. As a result, 
more shallow crack trajectories are predicted in the present study. Neverthe- 
less, the computed fracture parameters are comparable with previous results; 
thus, a similar fatigue life is anticipated. 

7.2.6 Concluding Remarks for Curvilinear Crack Growth 
Simulation in Narrow Body Fuselage Panel 

Curvilinear crack growth in a generic narrow body fuselage was studied. Com- 
parisons with experimental measurements suggest that the fracture toughness or- 
thotropy plays an important role in predicting the early crack growth trajectories. 
The subsequent crack growth after the initial crack deflection followed a trajectory 
where the local stress states are of a Mode I type. Thus, like crack growth in DCB 
specimens, crack turning and flapping as the crack approaches the tear strap is 
highly related to the T-stress. 



Stress Intensity Factor (ksi|\/inch) 
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Figure 7.15: Computed stress intensity factors versus half crack extension. The 
hollow and solid markers denote the computed stress intensity factors 
from the current isotropic prediction with r c = 0 and those from Chen 
et al. [21], respectively. 
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Figure 7.16: Computed stress intensity factors versus half crack extension. The 
hollow and solid markers denote the computed stress intensity fac- 
tors from the current isotropic prediction with r c = 0 and those from 
Potyondy et al. [106, 108], respectively. 
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The predicted results based on the proposed methodology show the potential 
to characterize curvilinear crack growth, but further studies as discussed in Sec- 
tion 7. 2. 4. 2 need to be conducted to fully assess its applicability as part of a damage 
tolerance methodology. 

7.3 Discussion: Elastic-Plastic Curvilinear Crack 
Growth and Residual Strength Prediction 

The above curvilinear crack growth simulations explicitly assume that the crack is 
grown under small scale yielding conditions. To simulate elastic-plastic curvilinear 
crack growth, one can use the predicted curvilinear crack path as the predefined 
crack path and stable crack growth and residual strength analyses can be performed 
accordingly. The procedure was used in Chen et al. [21, 20] to study the trajectory 
effect on residual strength prediction. 

A more rigorous approach is to use a directional criterion based on the current 
elastic-plastic states at crack tips directly. A procedure for mapping the state 
variables from one finite element mesh to another is then performed as the crack 
propagates. A plane stress, non-self-similar elastic-plastic crack growth simulation 
based on the CTOD directional criterion discussed in Section 5.4.2 in conjunction 
with the CTOA crack growth criterion has recently been implemented [62], The 
predicted results are comparable to those observed in the Arcan fracture tests. 
Future work is needed to assess the applicability of the mapping algorithm and 
direction criterion to fuselage structures under stable tearing. 

7.4 Summary 

The directional criteria developed in Chapter 5 were used to predict curvilinear 
crack growth in coupon tests and in full-scale fuselage panel tests. The predicted 
trajectories were in good agreement with those observed in the tests. 

The influence of various parameters on the crack trajectory prediction was 
studied. Both T-stress and fracture orthotropy were found to be essential to predict 
the observed paths. The proposed methodology shows its potential to predict crack 
turning and flapping that can be used as part of a damage tolerance methodology. 



Chapter 8 


Summary, Conclusions, and 
Recommendations for Future 
Work 


This chapter summarizes the contributions of this thesis, draws conclusions, and 
where appropriate, provides recommendations for future work. 

This dissertation begins with a description of aging aircraft problems faced 
by the aircraft community. Aging of aircraft may significantly reduce structural 
integrity and residual strength below an acceptable level. This concern serves 
as the primary motivation for the dissertation. The objective is to develop an 
accurate structural analysis methodology and a useful and usable software tool for 
predicting the structural integrity and residual strength of pressurized, thin-shell 
structures. 

Background material related to structural integrity of aircraft fuselages and 
effective simulations of arbitrary crack growth is discussed in Chapter 1. This 
serves as a departing point to study simulations of fracturing processes in thin- 
shell structures. The dissertation is then divided into two parts to facilitate the 
discussion. Chapters 2 through 4 deal with the crack tip opening angle (CTOA) 
fracture criterion obtained from coupon tests to the prediction of fracture behavior 
and residual strength of built-up aircraft fuselages. Chapters 5 through 7 discuss 
relevant issues for crack trajectory prediction methodologies to improve structural 
integrity of airframes. Summaries, conclusions, and recommendations for future 
work of each part are presented below. 
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8.1 Part One: Elastic-Plastic Crack Growth Sim- 
ulation and Residual Strength Prediction 

8.1.1 Chapter 2: Theory for CTOA-Driven Elastic-Plastic 
Crack Growth and Residual Strength Analysis 

Various fracture mechanics methods for simulating elastic-plastic crack growth and 
predicting residual strength of thin-sheet metallic structures are reviewed and cri- 
tiqued in Chapter 2. The fracture analysis methods include linear elastic fracture 
mechanics (LEFM) and elastic plastic fracture mechanics (EPFM) versions of K R , 
G r , J r , 5 r , T* r , and CTOA using two-dimensional and three-dimensional anal- 
yses. Among the methods, the CTOA fracture criterion with three-dimensional 
elastic-plastic analyses is found to be superior because of its relative independence 
of the geometry of the structure, the length of the crack, and the presence of mul- 
tiple cracks. Elastic-plastic crack growth, link-up of multiple cracks, and residual 
strength analyses using the CTOA fracture criterion are discussed. 

8.1.2 Chapter 3: Residual Strength Analysis of a Flat Panel 
with Self-Similar Elastic-Plastic Crack Growth 

In Chapter 3, numerical simulations of flat panel tests are conducted by using thin- 
shell finite element analyses. The CTOA fracture criterion is used to characterize 
elastic-plastic crack growth. Two sets of fracture tests are simulated: one with a 
single crack but different widths and the other with multiple cracks. 

Predicted results of the flat panel simulations with a single crack show two dis- 
tinct failure phenomena. For small specimens, plastic zones reach the free bound- 
ary and the limit load is attained due to net section yielding. For large speci- 
mens, plastic zones are well confined by the elastic region and residual strength is 
reached due to fracture instability. Results of predicted residual strength are com- 
parable to experimental measurements. Yet as the width of the panel increases, 
the relative difference between experimental measurements and numerical predic- 
tions increases. This discrepancy is associated with the three-dimensional nature 
of the stresses around the crack tip, a result of constraint effects due to the finite 
thickness of the panels. A plane strain core concept is proposed to incorporate the 
three-dimensional constraint effects into two-dimensional as well as thin shell anal- 
yses. Predicted results with the plane strain core follow those of three-dimensional 
analyses and experimental measurements for small and large panels. 

Predicted residual strength of small flat, panels with multiple cracks agrees well 
with experimental measurements. A loss of residual strength due to the presence 
of multiple small cracks is observed. 
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8.1.3 Chapter 4: Residual Strength Analysis of Fuselage 
Structures with Self-Similar Crack Growth 

Chapter 4 examines the feasibility and validity of the analysis methodology to 
predict residual strength of pressurized fuselage structures that are subjected to 
widespread fatigue damage (WFD). The first part of the chapter uses a relatively 
simple built-up configuration to examine the effect of lead crack sizes, multi-site 
damage (MSD), and material thinning due to corrosion damage. Predicted results 
indicate a 21.8 to 28.0% loss of residual strength due to the presence of small MSD. 
Coupling of MSD and corrosion damage leads to the most severe damage scenario. 

The second part of the chapter describes analyses of fuselage panels tested in a 
wide body, pressure test fixture. The objective is to validate the analysis method- 
ology by direct comparison of numerical and experimental results. A global-local 
hierarchical modeling strategy is used to analyze the panel tests. This modeling 
strategy allows one to obtain sufficient accuracy of computed values with reason- 
able computer resources. 

Predicted stress distributions in the vicinity of the lap joint are compared with 
strain gage readings. Major results from the strain gage comparison are: 

• For global and local models of about the same coarse mesh density, the pre- 
dicted results converge quickly and agree with experimental measurements. 

• Results with a much higher mesh density that is suitable for stable crack 
growth analysis disagree with the rest of the numerical predictions and ex- 
perimental measurements. The discrepancy is related to the idealized rep- 
resentation of the two-noded spring element for the rivet connection. The 
problem is effectively overcome by generating distributed connections be- 
tween the two-noded spring element and the surrounding shell elements. 

Elastic-plastic crack growth analyses using the CTOA fracture criterion are 
conducted. Numerical results for the case with and without MSD are compared 
to experimental observations. Two key factors are found to be crucial for accurate 
prediction of stable crack growth and residual strength of the wide body panel 
tests. One is to incorporate the residual plastic deformation left by the fatigue 
crack growth, and the other is to consider the failure of other structural elements 
during stable crack growth. The specific highlights are: 

• The residual plastic deformation or the plastic wake from fatigue crack 
growth has a strong effect on stable crack initiation but a mild effect on 
residual strength prediction. Neglecting plastic wake effect leads to a totally 
erroneous prediction of the early stable crack growth. 

• The breakage of the inner tear strap, categorized as possible failure of other 
structural elements during stable crack growth, is crucial to residual strength 
prediction. For all the analyses conducted, the occurrence of the broken tear 
strap reduces the predicted residual strength by 24% to 30%. 
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Both observed and predicted results of the wide body panel tests again show 
substantial reduction of residual strength due to the occurrence of MSD. 

8.1.4 Recommendations for Future Work 

For cracks in a pressurized fuselage, the out-of-plane deformation or bulging at 
the crack edges is an essential characteristic feature of the displacement fields. 
The current analysis methodology assumes that the same critical CTOA deter- 
mined from flat panel tests with guide plates suffices to characterize the fracturing 
processes. To fully justify this assumption, laboratory tests that generate the out- 
of-plane deformation during stable crack growth need to be conducted 1 . Numerical 
simulations of these laboratory tests using geometrically and materially nonlinear 
thin-shell and three-dimensional crack growth analyses will shed new fight on the 
possible invariance of the CTOA fracture criterion to crack bulging. 

Constraint effects due to finite thickness of the panels are currently incorpo- 
rated into thin-shell finite element analyses by an ad-hoc fashion, that is, using 
a plane strain core concept along the tearing crack path. Fully three-dimensional 
analyses or mixed thin-shell and three-dimensional analyses can automatically cap- 
ture the three-dimensional constraint effect and eliminate the need for the plane 
strain core. A mixed model consisting of thin shell and three-dimensional ideal- 
izations as illustrated in Figure 8.1 seems to be a very attractive approach. Using 
this approach, fracture behavior around the crack tip region can be described ac- 
curately by three-dimensional analyses while thin shell idealizations may apply to 
the remote regions where the through thickness effect can be ignored. 

The current model does not faithfully represent crack growth in the vicinity 
of rivets. Distributed connections may be adequate to represent load transfer 
through the rivets, but may not have sufficient accuracy to characterize fracturing 
processes. Further laboratory fracture tests and analyses of various lap-jointed 
configurations need to be conducted to quantify its effect on stable crack growth 
and residual strength prediction. 

The current analysis procedures of incorporating the residual plastic deforma- 
tion for stable tearing do not include the effect of crack face contact. This leads 
to a much higher crack-opening pressure in comparison with 2D plane stress re- 
sults and laboratory observations. Further study is needed to quantify its effect 
on residual strength prediction. 

Widespread fatigue damage (WFD) has two subsets: multi-site damage (MSD) 
and multi-element damage (MED). The effect of MSD on residual strength can be 
analyzed and accurately predicted by the current methodology. The MED is yet to 
be explored rigorously. Also, the similar scenario including static or dynamic failure 
of other structural elements during stable crack growth needs further investigation. 
A proper mechanism to initiate and propagate damage in other structural elements 

x The MT specimen with the width larger than 24 inch but without guide plates 
seem to be a plausible candidate. 
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Figure 8.1: Illustration of mixed modeling of thin shell and 3D solids. 


needs to be included in stable crack growth analyses. 

Material thinning due to corrosion damage that may occur in aging aircraft is 
modeled through a uniform thickness reduction over the skin at the lap joint. More 
detailed analyses are needed to assess its applicability in characterizing corrosion 
damage. 

Finally, the methodology developed herein is mainly for thin-sheet metallic 
structures. Its applicability to thick, heavily loaded structures (for example, wings) 
or to different materials (for example, composites) is yet to be determined. 


8.2 Part Two: Curvilinear Crack Growth Simu- 
lation 

8.2.1 Chapter 5: Theory for Curvilinear Crack Growth in 
Planar and Thin Shell Structures 

This chapter begins with the motivation for using the crack turning phenomenon to 
improve the structural integrity of fuselage structures. To predict a crack trajectory 
that is not known a priori , a criterion for predicting the crack propagation direction 
is required. 

The maximum tangential stress theory is used as a starting point to evaluate 
the direction of crack growth. Full stress and displacement fields in two-dimensions 
and asymptotic fields in thin plates subjected to bending are outlined. 

The crack growth direction criterion based on the two-dimensional, linear elas- 
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tic crack tip fields is assumed to be sufficient to handle thin shell problems under 
the LEFM conditions. Singular as well as non-singular constant stress fields are 
included in evaluating the tangential stress near the crack tip. A directional crite- 
rion based on the maximum tangential stress up to the order of the T-stress term 
is derived. The specific highlights include: 

• The predicted propagation angle is determined based on r c , A/, A'//, and 
and T. A predicted propagation angle diagram is presented using a non- 
dimensional parameter T where T = (8 T \/2Tir c ) / (3 Kj) . 

• Under pure Mode-I and positive T-stress conditions, the crack path instabil- 
ity will occur when r c > (9A'|)/(1287 tT 2 ). 

• Under general mixed-mode conditions, the criterion predicts a bigger prop- 
agation angle under positive T-stress and a smaller angle under negative 
T-stress. The criterion reduces to the Erdogan and Sih criterion when T = 
0 or r c = 0. 

The criterion is then extended to include the effect of fracture toughness or- 
thotropy. A simple elliptical function is used to characterize the anisotropic frac- 
ture resistance in different material orientations. The effect of fracture toughness 
orthotropy ratio and the material orientation angle on the predicted propagation 
angle is examined. 

The rest of the chapter examines possible extensions of the current crack growth 
direction criterion to handle geometrically and materially nonlinear problems. 

For elastic deformations of arbitrary magnitude, the extensions rely on finding 
the Lagrangian counterparts of conservative integrals, well-defined under elastic 
states with infinitesimal deformations, in the context of finite elastic deformations. 
The fracture parameters are then related to these conservative integrals. 

For elastic-plastic problems two directional criteria are examined: one based 
on the HRR fields and the other based on the crack tip opening displacement 
(CTOD) concept. The study concludes that neither the HRR type extension of 
the maximum tangential stress theory nor the CTOD based directional criterion is 
currently sufficient to fully characterize the direction of elastic-plastic crack growth. 

8.2.2 Chapter 6: Numerical Evaluation of T-Stress 

Numerical methods to obtain accurate T-stress for two-dimensional as well as thin 
shell problems are the main theme of Chapter 6. The specific highlights are: 

• Two types of path independent integrals for T-stress evaluations are pre- 
sented: one based on the Betti-Ravleigh reciprocal theorem and the other on 
Eshelby’s energy momentum tensor. The analytical and numerical equiva- 
lence between the two is found. 
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• A novel error estimator for T-stress is proposed to quantify and assess the 
accuracy of computed values. Specifically, it is found that the error of the 
computed T-stress is proportional to the ratio of stress intensity factor di- 
vided by the square root of the characteristic dimension of the integration 
domain where the path independent integral is evaluated. 

• Numerical accuracy in evaluating T-stress using the path independent inte- 
gral method is assessed by highly accurate two-dimensional p- and hp-version 
adaptive finite element analyses. 

8.2.3 Chapter 7: Curvilinear Crack Growth Simulations 

Chapter 7 analyzes curvilinear crack growth in double cantilever beam (DCB) 
specimens and in full-scale narrow body fuselage panel tests. The specific highlights 
for curvilinear crack growth in the DCB specimens are: 

• Observations in the fracture tests indicate that the crack tends to turn away 
from its initial crack tip under pure Mode I conditions. The crack growth 
directional criterion, incorporating the T-stress effect, is capable of capturing 
the essence of crack turning under such circumstance. The predicted results 
with r c = 0.09 inch best correlate the experimental data for stable tearing. 
The predicted results with r c = 0.06 inch best correlate the experimental 
data for fatigue crack growth. 

• The fracture toughness orthotropy is introduced to explain the difference 
between the measured crack paths in the L-T and T-L orientations. The 
orthotropic directional criterion shows a promising nature to correlate the 
experimental data. 

The specific highlights for curvilinear crack growth in the full-scale narrow body 
fuselage panel tests are: 

• T-stress has a very mild influence on the early crack trajectory prediction 
because of its small magnitude. But as the crack approaches the tear strap, 
T-stress increases and plays an important role in the crack turning prediction. 
For the case with r c = 0.09 inch, a sharp turning caused by T-stress is 
predicted as the crack approaches the tear strap. 

• The fracture toughness orthotropy has a strong effect on the early crack tra- 
jectory prediction. The predicted crack trajectory, with 10% higher fracture 
toughness in the T than in the L direction of propagation, agrees well with 
that from the experimental measurements, before the crack approaches the 
tear strap. 
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8.2.4 Recommendations for Future Work 

The current methodology assumes that the maximum tangential stress directional 
criterion developed under the two-dimensional, LEFM framework can be directly 
applied to thin shell problems. Further crack trajectory study of fracture coupon 
tests and full-scale fuselage panel tests is needed to fully justify the assumption. 

A rigorous elastic-plastic directional criterion for non-self-similar stable crack 
growth simulations is yet to be found. A procedure that maps the state variables 
from one finite element mesh to another as the crack propagates is yet to be 
implemented into the FRANC3D/STAGS software program. 

Accurate stress intensity factor and T-stress evaluations as the crack approaches 
the tear straps are crucial to predict the crack turning. In the current study, the 
convergence study was conducted to ensure the accuracy of fracture parameter 
evaluations. Further study on adaptive and higher order shell finite element anal- 
yses may help to improve the accuracy of numerical computation. Other numerical 
methods, for example, path independent integrals for geometrically nonlinear shells 
may also help to improve the accuracy of fracture parameter evaluations. 

The physical meaning of the parameter r c is yet to be found and the appropriate 
experimental method to measure r c is yet to be determined. Further understanding 
of fracture behavior at the meso- or micro-scale may shed new light on r c , and 
furthermore, the crack growth directional criterion. 
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